Factors Inﬂuencing Colonization and Survival of Juvenile Blue Crabs Callinectes sapidus in Southeastern U.S. Tidal Creeks

: Tidal creeks along the southeastern U.S. and Gulf of Mexico coastlines provide nursery habitats for commercially and ecologically important nekton, including juvenile blue crabs Callinectes sapidus , a valuable and heavily landed seafood species. Instream and watershed urbanization may inﬂuence the habitat value that tidal creeks provide to blue crabs. We investigated natural and anthropogenic factors inﬂuencing juvenile blue crab occupancy dynamics in eight ﬁrst-order tidal creeks in coastal North Carolina (USA). An auto-logistic hierarchical multi-season (dynamic) occupancy model with separate ecological and observation sub-models was ﬁtted to juvenile blue crab presence/absence data collected over replicate sampling visits in multiple seasons at three ﬁxed trapping sites in each creek. Colonization and survival are the processes operating on occupancy that are estimated with this formulation of the model. Covariates considered in the ecological sub-model included watershed imperviousness, the percent of salt marsh in each creek’s high tide area, percent salt marsh edge, site-level water depth, and site-level salinity. Temperature, salinity, and dissolved oxygen were covariates considered in the observation sub-model. In the ecological sub-model, watershed imperviousness was a meaningful negative covariate and site-level salinity was a positive covariate of survival probability. Imperviousness and salinity were each marginally meaningful on colonization probability. Water temperature was a positive covariate of detection probability in the observation sub-model. Mean estimated detection probability across all sites and seasons of the study was 0.186. The results suggest that development in tidal creek watersheds will impact occupancy dynamics of juvenile blue crabs. This places an emphasis on minimizing losses of natural land cover classes in tidal creek watersheds to reduce the negative impacts to populations of this important species. Future research should explore the relationship between imperviousness and salinity ﬂuctuations in tidal creeks to better understand how changing land cover inﬂuences water chemistry and ultimately the demographics of juvenile blue crabs.


Introduction
Development is increasing from 300-600% faster than human population growth in United States (U.S.) coastal zones [1], with these impacts impeding the ability of estuaries to serve as critical nurseries for resident and transient nekton species [2]. It is estimated that half of salt marsh habitats has been lost in the U.S. to urbanization of coastal seascapes [3]. Development in coastal areas also severs aquatic habitat connectivity [3][4][5] and reduces secondary biological production [5][6][7][8].
Tidal creeks along the southeastern U.S. and Gulf of Mexico coasts provide nursery habitats for the blue crab Callinectes sapidus. This species is a commercially and ecologically important crustacean found in the western Atlantic [21] that uses estuaries during its juvenile and adult phases (reviewed by [22]). The blue crab is the most important commercial species landed in North Carolina by ex-vessel value and by weight of its landings [23]. Salt marshes, which dominate the intertidal area of undeveloped tidal creeks in this region, are considered important secondary nursery habitats (i.e., second habitat used after initial settlement habitat) for juvenile blue crabs (>~20 mm carapace width (CW)) [24][25][26]. Small lower salinity estuaries may be particularly important nursery areas for this species in the southeastern U.S. [26].
The state of North Carolina has designated primary nursery areas (PNAs) to protect sensitive estuarine habitats from aquatic (e.g., trawling, dredging, shoreline armoring) and terrestrial impacts (e.g., elimination of upland buffers) in order to preserve their ecological value to fishes and crustaceans. PNAs include many first-order tidal creeks that are at risk of impacts from urbanization [27,28]. The impacts of urbanization to nekton inhabiting small tidal creeks either inside or outside of PNAs have not been extensively researched. One study in the southeastern U.S. region concluded that impacts to tidal creeks could be determined by a composite measure of urbanization-watershed imperviousness-and predicted negative effects to nekton will occur when imperiousness exceeds from 20-30% in tidal creek watersheds [15]. Using fishery-independent trawl data, Luczkovich et al. [29] found that increasing anthropogenic land use and changes in relative abundance of blue crabs in North Carolina PNAs were negatively correlated. However, to our knowledge, the effects of watershed imperviousness have not been examined using modeled demographics for any crustacean species. We explicitly model the demographics of blue crabs in this study to help reduce bias in estimated parameters and better understand how urbanization affects colonization and survival in tidal creek habitats.
There have been requests for more research to identify valuable habitats for juvenile crustaceans in North Carolina estuaries [27,30]. While there have been multiple studies evaluating estuarine habitat requirements of juvenile blue crabs [26,[31][32][33], they have typically evaluated habitat preferences from count data (e.g., catch per unit effort) analyzed outside the framework of hierarchical modeling. Indices of presence or abundance used to identify important habitat factors will bias estimates of demographics if detection probability is not accounted for in modeling [34,35]. In contrast, hierarchical models account for imperfect detection probability (false-negatives) typical of ecological count data, and thus, are capable of reducing the bias associated with the demographic processes that are modeled. One such family of hierarchical models is occupancy modeling, which estimates colonization and survival probability (or occupancy probability) between multiple seasons through model fitting to replicated presence/absence data on unmarked individuals while accounting for imperfect detection probability [36].
In this study, we fitted a multi-season (dynamic) occupancy model to data collected on juvenile blue crabs across a gradient of urbanization in North Carolina tidal creeks. Given increasing urbanization along the southeastern U.S. and Gulf of Mexico coastlines, this analysis helps determine how an important crustacean responds to urbanization and natural habitats within tidal creeks and their watersheds.

Study Sites and Field Sampling
Eight first-order tidal creeks were sampled in coastal North Carolina (USA) (Figure 1; Table 1). The studied creeks typify the patchiness and heterogeneity of estuarine habitats embedded within developing coastal landscapes [37] and possess a range of instream and watershed habitats and impacts representative of these systems along the southeastern U.S. [15] and Gulf of Mexico coastlines [7,38]. Except for their channels, undeveloped first-order tidal creeks in this region are largely intertidal, with this intertidal area dominated by the emergent vegetation Spartina alterniflora. One study creek can be considered undeveloped (Porters Creek: Figure 1, Panel F) in that there is no development in its intertidal or subtidal area except for a culvert at the upstream extent of tidal influence that is above all the sample sites in that creek. In contrast, developed first-order tidal creeks here are typically fringed with an armored shoreline that eliminates salt marsh [15,38,39]. Creeks sampled in this study were largely intertidal except for the completely sub-tidal and unvegetated Webb Creek (Figure 1, Panel B). The eight creeks are located approximately equidistant from the nearest ocean inlet through which ingress of blue crab megalopae occurs (Beaufort Inlet: 34.69 • N, 76.67 • W). Tidal creeks in this region lack submerged aquatic vegetation [26] but experience seasonal blooms of algae (e.g., Ulva lactuca) that may provide an ephemeral habitat for juvenile blue crabs [40]. Scattered clusters of oysters Crassostrea virginica are found in each creek. embedded within developing coastal landscapes [37] and possess a range of instream and watershed habitats and impacts representative of these systems along the southeastern U.S. [15] and Gulf of Mexico coastlines [7,38]. Except for their channels, undeveloped first-order tidal creeks in this region are largely intertidal, with this intertidal area dominated by the emergent vegetation Spartina alterniflora. One study creek can be considered undeveloped (Porters Creek: Figure 1, Panel F) in that there is no development in its intertidal or subtidal area except for a culvert at the upstream extent of tidal influence that is above all the sample sites in that creek. In contrast, developed first-order tidal creeks here are typically fringed with an armored shoreline that eliminates salt marsh [15,38,39]. Creeks sampled in this study were largely intertidal except for the completely sub-tidal and unvegetated Webb Creek ( Figure 1, Panel B). The eight creeks are located approximately equidistant from the nearest ocean inlet through which ingress of blue crab megalopae occurs (Beaufort Inlet: 34.69° N, 76.67° W). Tidal creeks in this region lack submerged aquatic vegetation [26] but experience seasonal blooms of algae (e.g., Ulva lactuca) that may provide an ephemeral habitat for juvenile blue crabs [40]. Scattered clusters of oysters Crassostrea virginica are found in each creek.   figure) and creek-specific maps (bottom of figure) for eight firstorder tidal creeks in coastal North Carolina (USA) studied for juvenile blue crab Callinectes sapidus occupancy dynamics. The medium and dark gray shading in each creek-specific map represent watershed area and impervious area within the watershed, respectively. The light gray shading and white in each creek-specific map denote high tide wetted area and permanently wetted area, respectively. The black dots in each panel represent fixed trapping sites (three per creek). Each of the three fixed sites per creek were sampled for juvenile blue crabs with a Gee-style 6-mm square-mesh minnow trap that was deployed over 33 replicate sampling events between summer 2018 and autumn 2020 ( Table 2). The only celestial season not sampled over this period was summer 2019. One sample site in each creek was located at or near the upstream extent of tidal influence while the other two sites were respectively located near the confluence of each creek with its higher-order estuary (Figure 1). Where a culvert was found, the two downstream sites were located either side of it (there were no culverts downstream of any of the study creeks). Mid-tide water depth at these sites ranged from 0.2 to 1.7 m and averaged 0.6 m (S.D. 0.3 m). Each trap was baited with frozen Atlantic menhaden Brevoortia tyrannus and fished sub-tidally over 24-h soak times near the sub-tidal/intertidal interface except for Webb Creek where traps were deployed in permanently sub-tidal habitats. For each replicate sampling event, all traps were set and then retrieved on the same respective days across all 24 sites. Trap catches were assumed to be independent among sites. Upon retrieval of each trap, blue crabs were counted, measured for carapace width (CW: mm), and released. Water temperature ( • C), salinity (parts-per-thousand (ppt)), and dissolved oxygen (mg/L) were collected in conjunction with each trap deployment using a YSI Pro 2030 instrument (YSI Inc., Yellow Springs, OH, USA).

Dynamic Occupancy Model to Describe Juvenile Blue Crab Habitat Use
Royle and Dorazio [41] describe a multi-season (dynamic) occupancy model fitted through Bayesian inference. Like other occupancy models, the dynamic version contains a hierarchy for the state (ecological) process and the detection (observation) process that is contingent on the latent state [42]. Thus, the observation sub-model estimates the probability of detection when a site is occupied. The use of colonization and survival (or occupancy), rather than proxies for abundance (e.g., catch-per-unit-effort), as state variables, is a recommended approach when describing population dynamics using count data that includes infrequent encounters of the studied species [43]. What separates occupancy and other hierarchical models from more traditional models of wildlife population dynamics is the explicit modeling of the observation (detection) process that accounts for falsenegatives [36,44,45]. Detection probability, if not accounted for in the modeling framework, becomes confounded with estimates of state parameters in count surveys of populationlevel processes [36,[45][46][47][48]. The dynamic occupancy model is akin to a robust design [49] in that it assumes closure at the site level among replicate surveys within primary sampling periods but relaxes this assumption between successive primary periods. Designating celestial season as a primary period was a logical choice for dynamic occupancy modeling in this study given that the blue crab displays seasonal pulses of megalopa settlement in estuaries and subsequent juvenile recruitment to nursery habitats (reviewed by [50]). After an initial Bernoulli trial of occupancy in the first study period, occupancy in later periods depends on whether a site was occupied previously and two components of occupancy: colonization and survival [45]. The ecological sub-model for occupancy probability ψ at site i and season t describes the latent occurrence state z i,t as: while the observation sub-model for detection probability p at site i, replicate survey j and season t describes the observation process y i,j,t conditioned on the latent state: Thus, an estimate of detection probability is made possible in occupancy models by replicated sampling within each primary period [42].
A temporal auto-logistic parameterization of the dynamic occupancy model was fitted to blue crab data in this study; this parameterization specifically models the colonization and survival processes that influence occupancy probability [41]. The parameter a in this logistic model is related to colonization, γ, by the equation: while the parameter b is related to survival, Φ, by the equation: The auto-logistic parametrization can accommodate environmental effects as modeled covariates. We fitted an auto-logistic parameterization of the dynamic occupancy model that included separate effects of each modeled covariate on colonization and survival. For example, a model that includes the effect of a covariate x on the probability of colonization or survival, takes the form: where π i,t represents colonization or survival for site i when transitioning to season t. A site that is previously occupied in season t − 1 can only have survival while a site that was not previously occupied can only have colonization. The regression coefficient α 1 represents the logit-scale influence of the covariate x on colonization and the sum of α 1 and α 2 represents the logit-scale influence of the covariate x on survival. Thus, each covariate in the ecological sub-model can be considered for its effect on colonization and survival.
Covariates can also be incorporated into a regression model describing the observation process in a hierarchical modeling framework. For example: describes a regression model where the logit-transformed detection probability p at site i, replicate visit j, and in season t is a function of the intercept β 0 and the regression coefficient β 1 of covariate x that was measured at site i, replicate visit j, and in season t. If measured during each repeat sampling event, the same environmental covariate (e.g., water temperature) can be considered in both the ecological and observation submodels. However, the influence of a covariate in the ecological sub-model is evaluated on a spatial scale (i.e., has one value per site) while it is evaluated on a spatio-temporal scale in the observation sub-model (i.e., has a unique value for each repeat sampling visit to a site).
Occupancy was modeled using the auto-logistic formulation described above given that our motivation was to study the influences of environmental effects as covariates of occupancy dynamics of juvenile blue crabs in tidal creeks. Covariates considered in the ecological sub-model included percent watershed imperviousness, the percent of each creek's high tide wetted area comprised of salt marsh, and the percent of linear vegetated edge relative to marsh area found in each creek. These landscape-scale factors (imperviousness, percent salt marsh, percent salt marsh edge) had equal values among the three sites in each creek. Details of how each of these variables were measured are found elsewhere [51]. Briefly, for watershed imperviousness, ArcGIS software (version 9.3.1) was used to estimate the watershed area of each creek and the percent of impervious surface within each watershed. Data used for this analysis were digital elevation models (DEMs) developed from light detection and ranging data, field survey data, and aerial photographs. Each watershed was estimated using the 'watershed' spatial analysis tool in the hydrology toolset in ArcToolbox ® . Extraction of relevant features (hill slopes and flow paths) from the DEMs was performed using the 'flow direction' tool, creating a new layer where flow direction and channeling points were represented. This layer was then used as the input for the 'watershed' tool which created a layer delineating watershed. From this layer, the watershed area could be estimated from the layer's attributes. Satellite imagery from Google Earth ® (images captured in 2010) was overlaid on watershed areas and used to estimate land use within the watershed. Impervious surface was estimated by manually drawing polygons over impervious surface areas within each watershed. Percent salt marsh and percent edge were measured upstream from the intersection of each creek with its higher-order estuary. Watershed imperviousness was modeled because it is considered a composite measure of anthropogenic impacts to tidal creeks in the southeastern U.S. [9,15] while salt marsh is the dominant habitat-type in this region's tidal creeks that lack instream impacts [10] and provides habitats for juvenile blue crabs [25,[52][53][54][55][56]. Percent edge was incorporated into models because it represents a metric of accessibility to the marsh surface [37,57]. Water depth was considered in the ecological sub-model because it structures nekton habitats in tidal creeks [58] and can influence mortality on juvenile crustaceans in shallow-water estuaries [32]. Whether each site was above or below a culvert was considered as a factor in the ecological sub-model owing to the potential for culverts to inhibit nekton movement in tidal creeks [19]. We did not include shoreline armoring in the ecological model; while shoreline armoring represents a common anthropogenic impact to nekton communities in southeastern U.S. tidal creeks [59], its presence at study sites was highly negatively correlated with the presence of salt marsh at these sites (Spearman rho = −0.80, p < 0.001). Taken together, these effects represent common natural and anthropogenic features within tidal creeks in the southeastern U.S. coastal region [15,39,51]. Finally, mean salinity was considered in the ecological sub-model because salinity can influence juvenile blue crabs in estuaries [26]. For salinities in the ecological sub-model, we used the mean of all salinity measurements taken at a particular site obtained from sampling throughout the study. We considered only the colonization effect of the culvert location in the ecological sub-model. For each of the other effects considered in the ecological sub-model, we considered both colonization and survival effects. We did not fit the ecological sub-model to water temperature given that changes in estuarine water temperature are inherent to changes in the celestial season at this latitude.
Given that blue crabs can respond to changes in estuarine water quality (reviewed by [22]), an observation sub-model was fitted with water temperature, salinity, and dissolved oxygen considered as regression covariates of detection probability. The observation sub-model was fitted to covariate data that were collected during each replicate sampling visit. In the case of passively fished traps, detection probability estimated via the observation sub-model represents the probability of capturing the target species if it is present in the vicinity of the sampling gear and estimated during periods of assumed closure (i.e., within each season).
The overall model was fitted to a three-dimensional array of observed presence/absence (binary) data where the number of sites, maximum number of replicates within a single season (seven), and number of consecutive seasons over the full duration of the study (10) represented each dimension. The array of response data was populated with 'NA's for the single unsampled season within the study as well as in instances when less than the maximum within-season number of replicate samples was collected; matrix data populated with 'NA's do not contribute to the likelihood. Each covariate in the ecological sub-model was fitted as one-dimensional vector consisting of a single value per site. Water temperature, salinity, and dissolved oxygen data were each incorporated into the observation sub-model as a three-dimensional array with equal dimensions as the array of trapping data. Each continuous covariate was standardized to increase the efficiency of model convergence [60].
Fitting used a customized version of published probabilistic code for the dynamic occupancy model [41,45] (Supplementary Material S1). An uninformative uniform prior probability distribution ('prior') with minimum and maximum values of 0 and 1, respectively, was assigned to the initial occupancy probability parameter ψ. An uninformative normal prior (mean and precision values of 0 and 0.01, respectively) was used for each of the two auto-logistic parameters and for each partial regression coefficient considered in the ecological and observation sub-models. For the culvert categorical factor, the absence level was set to zero and the effect of presence was evaluated relative to it ('effects parameterization': [61]). Initial values for each modeled parameter were generated over the range of its prior using the built-in rnorm function in R [62].
The model contained multiple calculated values (Supplementary Material S1). These included calculated effects of each covariate retained in the final model on logit-scale colonization and survival probability. We also calculated back-transformed colonization and survival probability as a function of each covariate retained in the final fitted model. These calculations were made over the range of standardized covariate values measured in this study for each covariate retained in the model while holding standardized values of each of the other retained covariates at their means; using this approach, other covariates drop out of each calculation (the standardized value of a mean is zero).
Models were fitted by calling OpenBUGS software (version 3.2: [63]) from R using the interface package R2OpenBUGS [64] and updating 50,000 times with every iteration saved after discarding the initial 40,000 updates as burn-in. Covariates whose 2.5% and 97.5% quantiles for both their partial regression coefficients (effect on logit-transformed colonization) and calculated values for effect on logit-transformed survival that most broadly overlapped zero were sequentially eliminated from subsequent model runs; the final model run included only covariates with posteriors of whose posterior masses did not widely overlap zero for the effect of the covariate on logit-transformed colonization or survival. Stationarity of each modeled parameter was determined by computing a Gelman-Rubin statistic (R) for it [65] and examining software-generated trace plots of retained updates to insure the adequate mixing of model chains.

Results
Across the eight creeks and 24 sites, there were 792 trap sets over 33 replicate sampling events; 5, 8, 12, and 8 replicate sampling events occurred in the spring, summer, fall, and winter, respectively. A total of 176 juvenile blue crabs were captured in 118 trap sets giving a 14.9% observed occupancy. Average carapace width (CW) was 31.7 mm (S.D. 9.6; range from 10 to 57). The observed proportional occupancy of traps was highest in summer (18.2%), followed by spring (17.5%), autumn (13.5%), and winter (13.5%). Mean (S.D.) water temperature, salinity, and dissolved oxygen across all samples and sites were 20.5 (7.4) • C, 16.0 (9.5) parts-per-thousand, and 5.4 (2.1) mg/L. Water temperature and dissolved oxygen generally overlapped among the three sites in each creek (Table 3). In contrast, salinity values were generally lowest at the upstream site in each creek (Table 3). Table 3. Mean (± S.D.) values for water temperature ( • C), salinity (parts-per-thousand), and dissolved oxygen (mg/L) at each of the three fixed sample sites in eight first-order tidal creeks in coastal North Carolina (USA) where juvenile blue crab Callinectes sapidus occupancy dynamics were studied. Sites in each creek are numbered from lowest to highest values (1-3) in a downstream-to-upstream orientation. Creeks are tabled in order (from west to east) that they appear on the overview map (see Figure 1).

Creek
Site Observed proportional occupancy varied widely within seasons (Figures 2 and 3). The magnitude of colonization and survival by season are represented by the intercepts; values for these intercepts (a for colonization and a + b for survival) ( Table 4) were generally consistent among the nine different intervals between successive seasons in that the 95% quantiles for each of these estimated values overlapped zero (Table 4). When averaged for the intervals between successive celestial seasons, the magnitude of colonization and survival were generally consistent among seasons (overlapping credible sets among seasons: Table 5).
The final dynamic occupancy model fitted to the crab data included two meaningful covariates in the ecological sub-model (Table 4): watershed imperviousness and mean site-specific salinity. The bulk of the posterior distribution was negative for the effect of watershed imperviousness on colonization while the bulk of the posterior was positive for the effect of mean site-level salinity on colonization (Table 4). For survival, almost all the posterior mass for the effect of watershed imperviousness was negative and the entire posterior mass for the effect of mean site-level salinity was positive (Table 5). Modeled colonization and survival probabilities decreased at higher percentages of watershed imperviousness ( Figure 4) and increased at higher salinities ( Figure 5).  The final fitted occupancy model included a single meaningful covariate in the observation sub-model. Water temperature was a positive covariate of detection probability ( Figure 6, Table 4). The masses of the posteriors for the other two effects considered in the observation sub-model (mean site-level salinity and dissolved oxygen) widely overlapped zero. At the mean occasion-specific water temperature, the back-transformed estimate of detection probability across all seasons and sites was 0.186. The value for the convergence statistic (R) was less than 1.1 for each modeled parameter in both the preliminary and final models; these statistics and trace plots indicated the adequate mixing of chains of retained updates of each modeled parameter.

Cybersecurity Awareness Campaign
As described in the previous paragraphs, an intensive awareness campaign through the IT departments of Greek Healthcare institutions was initiated, in December 2019, focusing on actions and precautions that each healthcare employee should undertake to protect the data they handle. A variety of communication means were employed, including:

1.
A certified GDPR training program provided by the Greek National Centre of Public Administration and Local Government to public servants.

2.
A flyer with important cybersecurity notes and indications which was distributed to all departments and clinics. In compliance with the Directive (EU) 2016/1148 of the European Parliament and of the Council of 6 July 2016 concerning measures for a high common level of network and information systems security across the Union, the flyer informed its readers that healthcare organizations have to comply with certain cybersecurity rules regarding their network and information systems. Consequently, the healthcare workforce was advised to:    [5] logit-scale colonization probability between seasons 5 and 6 −3.27 9.65 24.24 intcol [6] logit-scale colonization probability between seasons 6 and 7 −8.74 5.20 18.86 intcol [7] logit-scale colonization probability between seasons 7 and 8 −6.01 9. 19 24.01 intcol [8] logit-scale colonization probability between seasons 8 and 9 −6.63 7.64 22.56  [2] logit-scale survival probability between seasons 2 and 3 2.17 10. 53 28.48 intsurv [3] logit-scale survival probability between seasons 3 and 4 3.73 11.99 29.42 intsurv [4] logit-scale survival probability between seasons 4 and 5 −15.15 3.95 27.04 intsurv [5] logit-scale survival probability between seasons 5 and 6 −6.51 13.66 32.64 intsurv [6] logit-scale survival probability between seasons 6 and 7 3.24 9.42 26.79 intsurv [7] logit-scale survival probability between seasons 7 and 8 3.40 13.53 34.48 intsurv [8] logit-scale survival probability between seasons 8 and 9 4.49 14.02 31.46 intsurv [9] logit-scale survival probability between seasons 9 and 10 −11.19 −3.93 3 posterior mass for the effect of mean site-level salinity was positive (Table 5). Modeled colonization and survival probabilities decreased at higher percentages of watershed imperviousness ( Figure 4) and increased at higher salinities ( Figure 5).   The final fitted occupancy model included a single meaningful covariate in the observation sub-model. Water temperature was a positive covariate of detection probability ( Figure 6, Table 4). The masses of the posteriors for the other two effects considered in the observation sub-model (mean site-level salinity and dissolved oxygen) widely overlapped zero. At the mean occasion-specific water temperature, the back-transformed estimate of detection probability across all seasons and sites was 0.186. The value for the convergence statistic ( ) was less than 1.1 for each modeled parameter in both the preliminary and final

Discussion
This research investigated factors that could influence the occupancy dynamics of juvenile blue crabs (10 to 57 mm CW) inhabiting first-order tidal creeks in coastal North Carolina. These creeks and their watersheds typify the variable impacts due to urbanization in the southeastern U.S. and Gulf of Mexico coastal zones. The elimination of salt marsh and changes to watershed land cover types are typical consequences of urbanization in this region [15,51,66]. Watershed imperviousness had largely negative effects on colonization and survival probabilities of juvenile blue crabs while site-level salinity had largely and entirely positive effects on colonization and survival probability, respectively. The results highlight the need to conserve natural land cover classes so that tidal creeks can continue to provide productive nursery value for nekton relative to historical (pre- Figure 6. Backtransformed estimated detection probability (y-axis) vs. water temperature ( • C) (x-axis) from fitting a dynamic occupancy model to study juvenile blue crab Callinectes sapidus demographics in eight North Carolina (USA) tidal creeks. The range of plotted water temperatures spans the range of measurements taken during replicated sampling visits to fixed sampling sites. The median estimated detection probability is shown by the black line while 97.5 and 2.5 quantile estimates are shown by the upper and lower gray lines, respectively.

Discussion
This research investigated factors that could influence the occupancy dynamics of juvenile blue crabs (10 to 57 mm CW) inhabiting first-order tidal creeks in coastal North Carolina. These creeks and their watersheds typify the variable impacts due to urbanization in the southeastern U.S. and Gulf of Mexico coastal zones. The elimination of salt marsh and changes to watershed land cover types are typical consequences of urbanization in this region [15,51,66]. Watershed imperviousness had largely negative effects on colonization and survival probabilities of juvenile blue crabs while site-level salinity had largely and entirely positive effects on colonization and survival probability, respectively. The results highlight the need to conserve natural land cover classes so that tidal creeks can continue to provide productive nursery value for nekton relative to historical (pre-development) conditions [67]. Interpretations of our modeling results and their relationships to prior findings are discussed in greater detail below.
The effects of watershed imperviousness on estuarine nekton inhabiting tidal creeks appears specific to the species and metrics studied [15,59,68]. Two watersheds in this study had imperviousness exceeding the upper limit for the 20-30% purported threshold for negative impacts on tidal creek nekton [15]; thus, our results support the viewpoint that watershed imperviousness exceeding this threshold will negatively impact nekton inhabiting tidal creeks. It is not clear what specific mechanism by which greater imperviousness reduces rates of juvenile crab survival; one of these mechanisms may be that greater imperviousness reduces and/or leads to higher variability in salinities in tidal creeks (below). The impact of imperviousness on blue crab occupancy dynamics should be informative to coastal planners faced with balancing development with preserving tidal creek habitats [69] for juvenile blue crabs. We are unaware of how these rates of colonization and survival of juvenile blue crabs would change in higher-order estuaries where the effects of imperviousness may be less acute. Biotas are apt to show the strongest and most direct responses to urbanization in first-order tidal creeks because these systems are the most immediate estuarine repositories of land-based runoff [9,66,[70][71][72].
Juvenile blue crabs appear to possess wide salinity tolerances [73]. We found that colonization and survival probabilities increased with increasing site-level salinities; the trend was most pronounced for the effect of salinity on survival and is consistent with increasing survival with higher salinities in a study of tethered juvenile blue crabs in the Cape Fear River Estuary, North Carolina [26]. Several studies have found that the highest juvenile blue crab densities occur in polyhaline waters [55,56,[74][75][76], which is also consistent with our modeling on the importance of higher salinities. While we did not sample or test for an effect of watershed imperviousness on variations in site-level salinity, the finding of lower colonization and survival at lower salinities has interesting ramifications for how watershed land cover types may influence tidal creek salinity and, in turn, demographics of juvenile blue crabs in tidal creeks. It has been found that impervious land cover in tidal creek watersheds exacerbate salinity fluctuations due to runoff from hardened surfaces [7,8,72]. Quantitatively examining the interaction between imperviousness and the frequency and intensity of precipitation for their effects on salinity in tidal creeks would be a useful topic of future research. Our salinity measurements only occurred during trap sampling and were not collected at a fine enough temporal resolution to capture short fluctuations in precipitation and associated drops in salinity.
Salt marsh coverage, a composite indicator of ecological services provided within intertidal systems [12,13,31,77], did not influence the probability of colonization and survival to the same degree as watershed imperviousness and salinity, though the bulk of its posterior mass was positive for its effects on colonization and survival probability. Although empirical studies have concluded that juvenile blue crabs prefer structured habitats [25,33,54,76,78], the tidally limited accessibility of the salt marsh [79] or the difficulty for juvenile blue crabs to navigate it [80] could have masked the influences of salt marsh coverage, salt marsh presence, and salt marsh edge on occupancy dynamics. Juvenile blue crabs have been found over unstructured mud and sand [33,54,81] and in some cases appear to prefer the latter habitat [80]. It is believed that unstructured habitats represent essential secondary nursery habitats for larger juveniles (>25 mm CW) that have emigrated from structured habitats that provide refuge and are considered primary nursery areas for smaller conspecifics [50]. Unstructured habitats in tidal creeks include deeper sub-tidal channels. A depth effect was not found in this study even though others have concluded that juvenile blue crabs find refuge from predation in shallow (<0.3 m deep) waters relative to deeper areas (>0.7 m) [82]. One possibility is that a depth effect co-varies with other factors not measured here, such as food availability. Thus, depth could influence colonization or survival when modeled alongside other factors that were not tested in this study. Further, culvert location relative to each site did not meaningfully influence juvenile blue crab occupancy dynamics. The lack of a culvert effect does not rule out their potential effect on the settlement rates of planktonic blue crab megalopa above these structures, if they restrict or modify tidal flows [19].
Untested effects, not studied here, could influence occupancy dynamics of juvenile blue crabs. For example, lower abundances of juvenile blue crabs have been found along armored shorelines relative to fringing salt marsh [83][84][85], with increased predation risk cited as a mechanism for this trend [86]. However, the result here is consistent with conclusions that juvenile blue crabs are habitat generalists (reviewed by [22]). Juveniles 20-50 mm CW use heterogeneous structured and unstructured habitats as they grow out of sizes most vulnerable to predation and then appear to prioritize habitat usage based on prey densities [50].
The observation sub-model tested for an effect of several water quality characteristics on the detection probability of juvenile blue crabs. Water temperature was a positive covariate of detection probability. Thus, juvenile blue crabs had a higher probability of capture at occupied sites when the water was warmer. The effect of water temperature is not surprising, given that blue crabs, being ectotherms, are more active at warmer temperatures [87]. Unlike water temperature, salinity and dissolved oxygen did not affect detection probability. While juvenile blue crabs can sense and avoid hypoxia [88], oxygen values were hypoxic (<2 mg/L), a level below which blue crabs negatively respond [89], in only five percent of measurements in this study. The observation sub-model tests for an effect of covariates on the probability of detection (capture) at sites. This is different than the ecological sub-model, which tests for the influence of one or more effects on occupancy of a site. Thus, unlike the results from the ecological sub-model, which found that more saline sites were more likely to be occupied, results from the observation sub-model indicate that salinity does not influence the detection probability of juvenile blue crabs at tidal creek sites that are already occupied.

Caveats and Potential Sources of Error
Several issues may have influenced the analysis into meaningful covariates of juvenile blue crab occupancy in tidal creeks. As mentioned above, untested factors, such as the dispersal of ingressing larvae or post-larvae [50,90,91] or tidal creek morphology [78], may influence juvenile blue crab occupancy dynamics. The study creeks are roughly equidistant from the nearest ocean inlet where megalopa ingress occurs and megalopa settlement can occur year-round at this latitude {50]. This provides one explanation for the similar patterns of colonization among celestial seasons in tidal creeks considered sites of secondary dispersal of juvenile blue crabs [50].
Different rates of predation on juvenile blue crabs are possible among study areas or habitat types [74,[92][93][94]. Specific sources of mortality, such as density-dependent cannibalism [22,75] and predation from estuarine fish species such as Scianeops ocellatus and Cynoscion nebulosus [32,92], could vary among creeks and as a function of the habitats available. Cannibalism of juveniles can arise from different sizes of conspecifics or from consumption of molting crabs by hard-shell conspecifics [75]. The quantity of potential predators, suitable prey, interspecific competitors for food, disease and parasitism, each of which were not measured in this study, are additional variables that may have varied among sites and may influence juvenile blue crab habitat preferences, and thus, detection, colonization, and survival in tidal creeks [22,95].
Sites were sampled with a passive gear-type due to the challenges in sampling with active gears in small, low-energy tidal creeks; these challenges include the limited number of sites to deploy and retrieve active gear in a standardized fashion and issues with gear avoidance/capture efficiency when actively sampling structured habitats. Passive gears have issues of their own, such as variable responses affecting catch rates due to activity levels of the target species. These considerations were built into the sampling protocol by standardizing soak times across seasons and incorporating water quality covariates, such as temperature, in the observation sub-model.
Detection probability of blue crabs from sampling with passively fished traps includes two processes: capture efficiency and recovery efficiency [96]. Capture efficiency is the proportion of targeted individuals in an area that enter the sampling device while recovery efficiency is the proportion of animals surveyed relative to those that entered; crabs that show an indifference or aversion to the trap or larger individuals that are unable to enter traps would be examples of reduced capture efficiency while small crabs egressing through the trap mesh or entrances, or within-trap predation, are examples of reduced recovery efficiency. Juvenile blue crabs experienced ontogenetic changes in habitat use [26,97] and sizes either not captured or not recovered by traps in this study may have had different patterns of occupancy than found here. Capture efficiencies of >25 mm CW blue crabs in trawls have ranged from 1.4% (in seagrass) to 25% in soft bottom [22] and this brackets the 18% detection probability we estimated in this study using a passive gear. Additionally, our sampling and modeling approach allowed us to determine the effect of water temperature on detection probability. Thus, water temperature should be considered when interpreting or applying published estimates of blue crab capture efficiency.
Describing population demographics through fitting hierarchical models is a recommended approach given that these models account for imperfect detection [61] which, if ignored, can affect the ability to determine meaningful environmental covariates of population-level processes [98]. However, like other models, occupancy models have assumptions [99] that, if violated, may bias parameter estimates. Dynamic occupancy modeling assumes that (1) the probability of initial occupancy and subsequent rates (colonization, survival) are constant among sites, or differences among sites are modeled using covariates, (2) there is no unmodeled heterogeneity in detection probability, (3) there is independence among survey outcomes that collect presence/absence data, (4) population closure within primary sampling periods, and (5) no false positive detections occur [99]. We built covariates into the ecological and observation sub-models to account for heterogeneity in occupancy dynamics as well as detection. Primary sampling periods are often seasons where potential violation of the closure assumption is low. Low rates of emigration have been observed for juvenile blue crabs inhabiting tidal creeks [78] and we assumed here that immigration to and emigration from specific sites during the primary periods was nil or low. Blue crabs experience pronounced seasonal pulses of settlement and recruitment to estuaries [25,80,93]; this is another reason that we chose celestial season, as opposed to year, as a primary period. Violations of the closure assumption would lead to biased estimates of occupancy probability [100]. Future work could test for this assumption. Lack of independence among samples at nearby sites would likely lead to occupancy probability (or its components) being overestimated [101]. The modeling assumption of no false positives was likely satisfied, given the taxonomic experience of the data collectors (authors of this study). Despite having these assumptions, the dynamic occupancy model, although requiring sampling over a longer period, can provide estimates of detection probability, colonization, and survival without the need for more involved experiments that may have their own set of artifacts (e.g., tethering, depletion studies).

Conclusions
Our findings suggest that the blue crab, a nektonic species in the U.S. Atlantic and Gulf of Mexico estuaries, experiences reduced survival with increasing development in first-order tidal creek watersheds. Given that urbanization in the southeastern U.S. is projected to at least double by 2064 [102], watershed imperviousness will likely have an increasing impact on populations of juvenile blue crabs in tidal creeks. Tidal creeks are often included in the many PNAs designated by the state of North Carolina for protection from certain instream development and commercial fishing activities. Thus, the findings from this study have ramifications for the quality of habitat that PNAs can provide if their watersheds are unprotected from development. Imperviousness can be considered an additional anthropogenic impact to biological productivity above and beyond the coast-wide losses of intertidal salt marsh via urbanization [3]. The findings here should guide natural resource planners who are challenged with offsetting the impacts of coastal development with mitigation strategies that can help preserve the habitat services that tidal creeks provide to nekton.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/d13100491/s1, S1. OpenBUGS software code used to fit an auto-logistic parameterization 1 of the dynamic (multi-season) occupancy model through Bayesian inference to study habitat and urbanization covariates of juvenile blue crab Callinectes sapidus occupancy dynamics in eight North Carolina (USA) tidal creeks. The full model (all covariates) is presented here. Effects not considered in the final model run are commented out with a single pound sign (#). Comments for the modeler are denoted by multiple pound signs.