Multi-scale spatial dynamics of the Chesapeake Bay nettle, Chrysaora chesapeakei

. Untangling organisms' multi-scale spatial distributions is challenging due to their interactions with environments at multiple spatial and temporal scales. We deployed an Adaptive Resolution Imaging Sonar (ARIS) in a Chesapeake Bay sub-estuary to investigate multi-scale spatial distributions of the bay nettle ( Chrysaora chesapeakei ) in May – September of 2016 and 2017. Nettles were found to be dispersed in aggregations of multiple individuals. The average density of bay nettles in 2017 was higher than in 2016. Small aggregations (<5 m) were persistent in both years but only contributed <10% of total observed nettles. Large patches (>30 m) contributed ~40% of the total observed nettles. Large patches were more common in creek habitat where nettle density was higher. Nettle density was found to hit a peak value once in 2016, while there were two density peaks in 2017. Different aggregation patterns were observed during the second peak period in 2017 in which the number of large patches increased dramatically. Within the surveyed waterscape, the spatial patterns were consistent over time with higher abundance in the source creek than in the river channel, which underscores that C. chesapeakei requires hard substrate in shallow creeks for its benthic polyp stage. Using the ratio between nettles in the creek and near creek mouth as a proxy for dispersal rate, more nettles were transported out of the creek in 2017 than in 2016. The increase in patch size and high dispersion rate in peak periods in 2017 suggests that individuals were moving out of the creek habitat as density increased. Results highlight the complex spatial structure of bay nettles, which has major impacts on density estimates and subsequently affects our understanding of jelly ﬁ sh population dynamics and long-term trends.


INTRODUCTION
Understanding spatial patterns and related ecological processes is a fundamental problem in ecology (Turner et al. 1989, Wiens 1989, Levin 1992, Rahbek 2005. However, quantifying spatial patterns at different spatial scales remains difficult. First, it is essential to sample the proper spatial extent because it is not feasible to infer spatial information out of sampling range. Second, proper sampling requires adequate resolution which is often determined by the life history and migration patterns of the target species (Henderson-Sellers et al. 1985, Meentemeyer andBox 1987). Third, spatially explicit processes at different scales may affect organisms' distribution at multiple spatial scales. For example, prey and habitat dynamics could have impact on predator distribution at multiple scales ranging from small to large scale (Benoit-Bird andAu 2003, Resetarits andSilberbush 2016).
Estuaries like the Chesapeake Bay exemplify the complexity of habitat structure and physical, chemical, and biological processes at multiple spatial scales (Costello et al. 2006, Vinagre et al. 2017. Undersampling is a common problem in shallow estuaries because, in general, traditional sampling methods like nets cannot adequately sample mobile and patchily distributed organisms due to gear avoidance (e.g., Gale and Mohr 1978) or clogging (e.g., Palmieri et al. 2014). Research vessels may have difficulty deploying gear or navigating narrow and shallow-water sites, and the complex habitat structures common in these areas are intrinsically difficult to sample (Rozas and Minello 1997, Bayley and Herendeen 2000, Breitburg and Burrell 2014. Surveys in many estuarine and coastal habitats, let alone shallow areas in these habitats, are infrequently optimized for medusae or other jellyfish life stages. Instead, jellyfish species are sampled as fortuitous additions in surveys designed to target fish or plankton (Boero et al. 2008, Duarte et al. 2013, Gibbons and Richardson 2013, Brodeur et al. 2016. Jellyfish predation affects zooplankton and ichthyoplankton, potentially with detrimental effects on fisheries. Partially sampling jellyfish habitat at a coarse resolution limits our understanding of the spatial dynamics of gelatinous species, resulting in an underestimation of their impact on ecosystem function and structure (Purcell andGrover 1990, Baird et al. 2004).
Spatial distribution and abundance of jellyfish are closely related to their life histories. The typical scyphozoan life cycle contains a free-swimming or pelagic medusa stage in addition to a sessile (benthic) polyp stage which is responsible for overwintering success and longevity in shallow habitat (Arai 1997, 2009, Albert 2005, Hamner and Dawson 2009, Lucas et al. 2012, Shahrestani and Bi 2018. Benthic polyps show preference for substrate type and many species readily proliferate on unnatural hard substrate, suggesting manmade structures such as pilings and hardened shorelines may promote polyp settlement (Brewer 1984, Pitt 2000, Holst and Jarms 2007, Hoover and Purcell 2009, Duarte et al. 2013). This could potentially lead to shifts in the abundance and biogeographic distribution of jellyfish species. Additionally, natural habitat loss caused by the destruction of oyster reefs may also affect the population and spatial dynamics of jellyfish species, as the polyp stage prefers hard substrate for settlement (Breitburg and Fulford 2006, Breitburg et al. 2010, Hubot et al. 2017. Within estuarine systems, settlement and growth of polyp colonies often occur in semisheltered, shallow areas (Cargo and Schultz 1966, Purcell and Grover 1990, Purcell 1992, Olesen et al. 1994, Breitburg and Burrell 2014, Shahrestani and Bi 2018. The dispersal of medusae from the creeks or semi-closed coves to the river channel or mainstem is facilitated by downstream tidal or wind-driven flow (Hartwick 1991, Moriarty et al. 2012.
Our study explored the localized patterns of Chrysaora chesapeakei distribution defined by a source-sink framework, highlighting the role of a small creek within the extent of a sub-estuary study area or waterscape (Patuxent River, Chesapeake Bay). To overcome the limitations of traditional sampling methodology, we deployed an Adaptive Resolution Imaging Sonar (ARIS) system to collect sonar data of C. chesapeakei in a Chesapeake Bay tributary and analyzed medusa density and distribution through space and time at a high sampling resolution. Our robust data collection and processing procedures allowed us to investigate spatial patterns at multiple scales with an unprecedented look at fine-scale C. chesapeakei spatial dynamics. We hypothesized that differences in the spatial dynamics of the Chesapeake Bay jellyfish reflect the interannual variation in abundance of C. chesapeakei noted between the 2016 and 2017 summer seasons. Quantifying and comparing variability in abundance and spatial patterns of medusae are important for understanding dispersal and subsequent life history success of C. chesapeakei (and other jellyfish spp.) in a dynamic shallow habitat where environmental conditions can vary greatly from year to year (Cargo and King 1990, Purcell et al. 2000, Breitburg and Fulford 2006.

Study sites
The Chesapeake Bay is the largest estuary in the United States, and its environmental conditions vary seasonally and interannually. The hydrodynamics of estuaries are tightly coupled with biogeochemical properties such as salinity, nutrient transport, and hypoxia (Hagy et al. 2000, Hong andShen 2012). We focused our study in the lower portion of the Patuxent River, a mid-bay tributary. We used discharge rates from the USGS Station 01594440 (https://waterda ta.usgs.gov/usa/nwis/uv?01594440) as a proxy to describe differences in winter and summer flow conditions of our Patuxent River waterscape ( Fig. 1) in 2016 and 2017 (Hagy et al. 2000). Wind direction and wind speed (mph) across the summers of 2016 and 2017 were sourced from Jefferson Patterson Park Station in St. Leonard, MD (www.wunderground.com).
Three transects were designed to reflect C. chesapeakei source-sink dynamics as they unfolded from source to dispersal habitat. The vertical St. Leonard's Creek (SLC) transect was approximately 3 km long, starting midway up the creek and ending at the mouth of the creek where it met the Patuxent River channel (Fig. 1). We continued the transect across the upper Patuxent River channel with a horizontal transect approximately 2.5 km in length (UPX). Our final transect was a horizontal cross section of the lower portion of the Patuxent River channel (LPX) and covered an average distance of 2 km. Within the sub-estuary and across our transects, we assumed jellyfish found in the Patuxent River channel were mainly sourced from St. Leonard's Creek (SLC) while the Patuxent River channel (characterized by UPX and LPX transects) was considered dispersal habitat. C. chesapeakei found in the lower Patuxent could also have been sourced from adjacent creeks including Cuckold Creek and Hellen Creek (Fig. 1).

Field sampling
Transects were sampled approximately biweekly from May 26 to August 28 in 2016 (8 sample events) and from June 08 to August 31 in 2017 (nine sample events) for a total of 17 cruises. Most sampling events began in LPX during high tide; then, we worked our way up to UPX and into SLC. We collected jellyfish spatial distribution data using an ARIS (Sound Metrics, Bellevue, Washington, USA) sonar imaging camera, which is capable of recording high-definition images in the characteristically high-turbidity waters of our study sites. The ARIS was attached to a polearm that pivoted on the gunwale of a~6 m research vessel (Fig. 2). A live feed was recorded using a laptop computer and ARIScope software (Sound Metrics) with an applied tow filter. A portable 125V generator powered both the sonar camera and the laptop computer controlling ARIS operations.
Once submerged approximately 1 m below the surface, the camera was towed at the slowest possible speed of 1 nautical mile/h, and data were recorded at a rate of 15 frames/s for the highest possible resolution (see example image in Video S1). In 2016, we deployed the ARIS at an average pitch of approximately 50°f rom the water's surface and pointed the camera toward the stern. In 2017, we increased the pitch to 75°and reversed the face of the camera to point ahead to the bow to improve image quality. In the first four sampling events of 2016, we subsampled the deep area of the river channel (UPX and LPX transects) by maximizing the camera imaging range to 15 m in hopes of detecting jellyfish at deeper ranges. Limited medusae were found at these depths (<5 medusae). Maximizing the imaging range limited resolution of the camera and made processing more challenging, so we opted to limit our sampling range to 7 m considering the low counts of medusae observed beyond this range. These changes were made to improve image quality and maximize the sampling volume. Temperature and salinity were recorded every second along the tow using an RBR Concerto CTD (RBR, Ottawa, Ontario, USA), which was submerged 1 m below the surface of the water. A handheld Garmin GPS unit (Garmin Fig. 2. The ARIS sonar camera (Sound Metrics) was towed along transect lines while attached to a small research vessel (A). The mounting gear included a pivot that attached to the gunwale of the research vessel and was fitted with a polearm that extends approximately 30cm into the water column. The ARIS camera was affixed to the end of the polearm with a mounting plate, and the pitch of the camera had an average deployment angle of 55°below the water's surface. Sound Metrics software (ARIScope) records footage as seen in panel B, where Chrysaora chesapeakei >5 cm in bell diameter were detectable and marked in the water column. The range markers denote the range of the field of view of the camera, and the dashed red line (B) is the range of an individual nettle as detected by ARISFish software. The surface area of the blue triangle (B) created by the compression of 3D data (ARIScope) was estimated as sample volume.
International, Olathe, Kansas, USA) was used to map the trajectory of the transects for the duration of each sampling event (~5 h).

Data processing
In 2016, we recorded more than 11.5 h of data (627,518 frames), and in 2017, we recorded 12.5 h Fig. 3. Conceptual diagram of transect binning, whereby data recorded along the transects were divided into 5 m. Within-bin density was estimated by counting medusae within a 5-m bin and estimating the volume of each 5-m bin. The focal length of the ARIS is 0.7-1.0 m, and organisms present in this volume of water are unsampled. In this image, these undetectable organisms are represented by red jellyfish. of data (676,113 frames) for a total of 1,303,631 frames. Each frame was manually inspected using ARISFish software, and positively identified jellyfish were counted and marked in the water column. Recorded depths of an individual were accurate to 0.1 m. Using interpolated transect tracks (s À1 ), we estimated the latitude and longitude of all counted C. chesapeakei. Sample volume and sample depth were calculated for each recorded frame using image analysis methods adapted from Shahrestani et al. (2017). Frames were aggregated to standardize jellyfish counts to metric units, whereby C. chesapeakei medusa occurrence was summed within cubic meter bins (Fig. 3). All counts were standardized by sample volume to account for variation in sampling effort. Count data are provided in Data S1.

Multi-scale spatial analyses
To investigate C. chesapeakei spatial patterns at 5-400 m, patches were enumerated for each cruise with at least 100 individuals observed in the surveyed area. A patch was a collection of adjacent 5-m bins along a transect with observed densities above the mean density for the transect and cruise (Fig. 4), assuming a Poisson distribution for the abundance. Patch sizes were estimated as the number of 5-m bins sampled within each patch and were estimated using the rle function in R.
To explore C. chesapeakei spatial patterns at 400 m-1 km along the distance of sampled transects, it was necessary to account for small-scale variability in the start and stop locations between cruises and among transects, as well as the meandering tracks in SLC due to the geomorphology of our study site (Fig. 1). To do so, we defined a variable "river distance" which represents the closest distance (m) an individual C. chesapeakei medusa is to a source point within a transect. The source point of each transect was defined as the northernmost (SLC) or easternmost (UPX, LPX) sample location on any given cruise within each of our three study sites. We investigated mesoscale patterns of spatial autocorrelation in the count data by binning the jellyfish observations into 5-m bins, using the river distance variable, and plotting autocorrelation functions (ACFs) for the data collected in all three transects in the peak week of each year.
Spatial difference among different habitats with the waterscape was quantified using a generalized additive model (GAM), which is a nonparametric extension of the generalized linear model. We chose to use a GAM as the modeling framework to incorporate the spatial and temporal autocorrelation which was inherent in the dataset. Model selection for the most parsimonious model included screening for low values of Akaike information criterion (AIC) as well as high values of global deviance (%) explained by the model. The final GAM was used to generate spatially explicit estimations of C. chesapeakei abundance within our study habitat across years.
In the spatially explicit GAM, we considered the spatial-temporal nature of our data by incorporating both a river distance variable, which accounted for spatial autocorrelation, and a week of the year variable, which modeled seasonality. The GAM was fit for both the creek (SLC) and channel (UPX and LPX) of the Patuxent River waterscape using nonparametric smoother functions with R statistical software (mgcv package). The following model was selected, where s represents a smoothing factor: Chrysaora chesapeakei count~negative binomial = year + s(week of year, by = creek/channel) + s(river distance, by = creek/channel) + log (volume).
For spatially explicit estimations, we used the raster package in R to grid and project our waterscape into 8 9 6.25 m grid cells, assigning both values of river distance and volume to each cell. The volume of each cell (400 m 3 ) within the waterscape was calculated by multiplying the surface area of each cell (50 m 2 ) by 8 m depth, whereby we limited predictions to the first 8 m of the water column due to the low number of medusae counted below that depth.

High-resolution sonar data and seasonal dynamics
We counted 57,846 jellyfish over the course of our study duration and mapped their geo-referenced locations in the water column to visualize their spatial distribution in the vertical and horizontal direction (depth vs. distance along each transect). The average density of C. chesapeakei was higher in 2017 (0.05 AE 0.22 individuals/m 3 ) than in 2016 (0.01 AE 0.10 individuals/m 3 ). The large standard deviation associated with the mean estimates suggested data were over-dispersed as indicated by the continuous counts along transects which revealed variable density distributions and dispersal patterns between sites and sampling events (Fig. 4). Average depth of observed medusae across sites, cruises, and years was between 2 and 4 m, although the average depth of C. chesapeakei was closer to the surface in 2017 (l = 2.3 m, SD = 0.8 m) compared to 2016 (l = 3.1 m, SD = 1.1 m).
The average density in SLC was higher than in the river channel (UPX and LPX) in both years (Table 1). However, the density difference between SLC and the river channel was much larger in 2016 than in 2017,~2 orders of magnitude higher in SLC during the July peak period. C. chesapeakei were more dispersed in the river channel than in SLC in both years as indicated by the coefficient of variation. However, the difference in dispersion between SLC and the river channel was larger in 2016 than in 2017 (Table 1).
The average density of C. chesapeakei showed large seasonal variation in both years, but with distinct different seasonal patterns (Fig. 5). In 2016, there was a single density peak on July 21 in the study area. In 2017, two density peaks were observed in SLC on July 18 and August 22.

Aggregations at 5-400 m
Across all transects, highest counts of medusae were found in the creek. In the source creek (SLC), 75% of the population was contained in <1.0% of the total volume sampled within each cruise, except 28 June 2016. During all peak density weeks, the densest aggregations occurred in SLC, where 75% of C. chesapeakei were contained in 0.6% of sampled volume in 2016 and in 0.2% of the sampled volume in 2017 (18 July, Fig. 4).
Regarding 2016 cruises after 28 June, 75% of the population was concentrated at the start of our SLC transect toward the headwaters (Fig. 1). This was not the case in 2017, where in four of the seven cruises the densest aggregations of C. chesapeakei were found dispersed~1.5 km downstream, including during the peak week of 18 July 2017 (Fig. 4).
Aggregations at 5-m scale were persistent throughout the entire study area, and the number of patches declined sharply as patch size increased (Figs. 6 and 7). Patch sizes in the SLC transect were often <50 m when C. chesapeakei was less abundant (Fig. 6). When C. chesapeakei became more abundant, particularly during peak periods, patch sizes ranging from 50 to 200 m were also common. In SLC, the largest patch size was observed during the July 2017 peak period. On 22 August 2017, patches <50 m were much less common than on other sampling dates and patches >50 m occurred more than other periods instead. In the Patuxent River channel, patch size typically ranged from 5 to 30 m when nettles were less abundant and could extend to 50 m during the peak period (Fig. 7). Patches at >200m scale did occur but were less common than in SLC.
The majority of observed C. chesapeakei occurred within patches. In the river channel, the number of C. chesapeakei within patches totaled 7218 individuals and made up 83% of total observed individuals in the river channel. In SLC, 25,457 individuals were observed within patches and made up 74% of total observed individuals in SLC. Although 5-m scale aggregations were persistent, the number of individuals only made up 10% and 5% of total observed individuals in the river channel and SLC, respectively. Large patches contributed more to the total number of observed individuals. In the river channel, individuals within patches >5 and >30 m made up 73% and 40% of the total observed individuals, respectively. In SLC, patches >5 and >50 m made up 69% and 45% of the total observed individuals, respectively.
The histogram of the total count in each 5-m bin along each transect in each year also suggested different aggregation patterns. The number of 5-m bins with >40 C. chesapeakei was highest on 18 July 2017 in SLC followed by the cruise on 21 July 2016 (Fig. 8). However, this was not the case for the second peak period of 2017 (around 22 August 2017), in which the number of bins with >40 nettles was relatively low and suggested a different aggregation pattern between years.

Spatial patterns at 400 m-1 km
The ACF plots showed strong patterns of sinusoidal decay, that is, density gradient of C. chesapeakei over space within the creek (Figs. 9A-C). However, the contrast between the peak abundance weeks of 2016 and 2017 was distinct. In 2016, the range of positive spatial autocorrelation was limited to 400 m and negative spatial autocorrelation started to occur as spatial distance increased to >400 m in SLC, which suggested that C. chesapeakei was more concentrated in upper than in lower SLC in 2016 (Fig. 4). In 2017, the spatial range of positive autocorrelation in both peak periods extended to 500-600 m which reflected the increase in patch size and density (Figs. 9B,C). The negative spatial autocorrelation after 500 -600m reflected the difference between upper SLC and lower SLC. The spatial patterns in the river channel showed distinct difference between 2016 and 2017. No distinct patterns were observed in the river channel in 2016 (Fig. 9D,G), but in 2017, clear spatial patterns were observed for both transects (Fig. 9E, H). The range of positive autocorrelation along the UPX transect was relatively short,~100 m, followed by negative spatial autocorrelation, which was likely caused by the plume at the SLC mouth (Fig. 4).

Habitat-related spatial patterns and dispersal
Water column depth, temperature, and salinity did not significantly explain patterns observed in the data as evidenced by performance of the selected model (GAM). The final model performed well and explained 84.8% (R 2 = 0.76) of the global deviance observed in the dataset. The final GAM had the lowest AIC value relative to other models which included salinity, temperature, and depth as possible explanatory values. A negative binomial distribution was selected as the best family for our data and greatly outperformed both Gaussian and zero-inflated Poisson distributions in the GAM. The spatially explicit model allowed us to estimate total abundance and investigate patterns of dispersion in our study habitat in the summers of 2016 and 2017 (Table 2). Estimated abundance of C. chesapeakei in the entire waterscape was approximately 39 higher in 2017 than in 2016 (Table 2). Higher density in 2017 than in 2016 was also apparent during peak abundance in the creek with 38 medusae/100 m 3 and 20 medusae/ 100 m 3 , respectively. The highest abundances predicted throughout entire waterscape occurred during the week of 21 July 2016 and 18 July 2017 at 495,179 and 1,609,360 individuals/m 3 , respectively. When investigating weeks with similar Fig. 7. Histogram of Chrysaora chesapeakei patch size in the Patuxent River channel for each cruise. Transects were binned into 5 m. A patch was a collection of adjacent 5-m bins along a transect with observed densities above the mean density for the transect and cruise, assuming a Poisson distribution for the abundance. abundance between years, we noted a stark contrast in the dispersion of jellyfish between our study sites. For example, 83% of the total observed jellyfish were in the creek in 2016, while in 2017, 88% of medusae were in the channel. Overall, abundance was higher in the creeks for 90% of the dates that were included in our computations, though abundance predicted in the channel during 2017 was higher than abundance predicted in the creek during 2016.

DISCUSSION
Disregarding jellyfish bloom impact on ecosystem function and production is a concern considering their high rates of consumption of zooplankton (competition) and ichthyoplankton (predation), which is thought to have had detrimental effects on fish populations worldwide (Cowan et al. 1992, Purcell et al. 1994, Suchman and Sullivan 2000, Purcell and Arai 2001, Sommer et al. 2002, Brodeur et al. 2008. In Chesapeake Bay, Chrysaora chesapeakei or the bay sea nettle is present during summer to early fall and contributes substantially to controlling populations of zooplankton (Purcell 1992) as a top predator (Clifford and Cargo 1978, Feigenbaum and Kelly 1984, Baird and Ulanowicz 1989, Sullivan and Kremer 2011. Results from our study revealed similar patterns of C. chesapeakei abundance and seasonality compared to previous studies conducted in the Patuxent River tributary (Cargo and Schultz 1966, Cones and Haven 1969, Loeb 1972, Calder 1974, Baird and Ulanowicz 1989, Purcell 1992, Suchman and Sullivan 1998, Brown et al. 2002, Breitburg and Fulford 2006, Decker et al. 2007, Breitburg and Burrell 2014, Tay and Hood 2017. Fig. 8. Relative Chrysaora chesapeakei density noted between years and across transects. Transects were binned into 5 m, and the relative density within each bin was categorized as low, medium, or high within years and each transect study site: St. Leonard's Creek (SLC), the upper Patuxent River channel (UPX), and the lower Patuxent River channel (LPX). In 2016 (A-C), low counts represent bins with 1-40 medusae, medium-density bins have 40-80 medusae, and high-density bins contain 81-120 medusae. In 2017 (D-F), the ranges are higher due to the higher overall density. Low counts represent bins that contain 1-60 medusae, medium bins have 60-120 medusae, and finally, high counts contain 120-180 medusae. Overall, the highest frequency of high-density nettle clusters occurred on 18 July 2017 in SLC. However, the difference in seasonal patterns between SLC and the river channels in 2017 was unexpected. In 2017, there were two peaks observed in July and August, in contrast to the single peak in 2016. The abundance was much higher in the SLC (source) than in the river channel (sink). Therefore, it is important to understand the impact of jellyfish blooms in the source habitat, where abundance tends to be much higher.

Aggregations and potential causes
Application of advanced technology has allowed for an unprecedented look at the spatial patterns of a jellyfish species in its vertical and horizontal distribution. In sampling thousands of jellyfish, we were able to identify the average depth of jellyfish in the vertical water column. In 2016 and 2017, we observed C. chesapeakei concentrations at midwater depths (2-8 m). Pycnocline aggregation is common with many scyphozoans species, where they are exposed to less turbulence, slower mass transfer, and higher prey concentrations (Keister et al. 2000, Graham et al. 2001, Purcell et al. 2014. While vertical distributions of jellyfish did not show stark variability between years, C. chesapeakei spatial patterns and abundance within and among our study sites were highly variable in 2016 and 2017 and much less understood.
Our results suggested <5 m aggregations were common throughout the season and the entire study, which may be explained by reproductive needs. The sexual production of planula larvae requires close proximity, as both males and females are broadcast spawners requiring proximity spawning at fine spatial scales (Arai 1997, Graham et al. 2001, Hamner and Dawson 2009 Aglieri et al. 2014). However, in some cases aggregations at scales <5 m could be less pronounced, such as in the second peak period of 22 August 2017. The causes for the different spatial patterns observed during the two peak periods were not clear, but physical forces such as windinduced wave and turbulence could be important factors (Arai 1992, Graham et al. 2001. Scale-variant spatial patterns are often detectable with species, such as regularity at small scales and clustering at larger scales (Jordano et al. 2003, Dale andFortin 2014). However, this phenomenon is difficult to capture with marine organisms, particularly in shallow-water habitat where sampling gear deployment and data resolution are limited. Using high spatial resolution data, we were able to detect spatial patterns at different scales over time. The spatial patterns at >50 m are likely a combination of high population density and potential physical drivers. The increase in patch sizes during peak periods, particularly the August peak period in 2017, may suggest density-dependent dispersal.
The occurrence of high-density aggregations of C. chesapeakei in the source creeks is likely correlated with a strobilation event (mass asexual reproduction) occurring 40-60 d prior, a time duration that allows for juveniles (1-2 mm) to grow to a size detectable by the sonar camera (the smallest bell diameter detected by the sonar camera was approximately 7-10 cm). In 2016, we assumed the appearance of medusae in the early summer in the creek reflected the first strobilation event, which would have occurred 40-60 d prior. This was followed by the appearance of the high-density aggregations, which could reflect a second strobilation event. From these observations, it would be reasonable to conclude there were only two possible strobilation events in 2016. However, in 2017, if we assume the appearance of C. chesapeakei as the first strobilation event and each of the two high-density peaks afterwards indicated one separate strobilation event, there were three strobilation events. This may begin to explain the high variability in bloom magnitude noted from 2016 to 2017. Strobilation events are the primary method of asexual reproduction and numerical growth for nonholoplankton scyphozoans worldwide (Arai 2009, Lucas et al. 2012, and the success of early life stages affects the magnitude of a jellyfish bloom from year to year (Cargo and King 1990). Though it is difficult to provide evidence for the exact mechanisms that lead to third strobilation event, the stark contrast in predicted abundance between years allowed for a unique opportunity to explore the multi-scale spatial patterns of C. chesapeakei medusae that appear to reflect density-dependent dispersal.
Spatial patterns and potential causes within the waterscape The settlement and perennial success of C. chesapeakei polyps in the Chesapeake Bay correlate significantly with residence time and salinity regimes, and sampling all regions important to annual processes is vital to obtaining the whole picture of nettle abundance and spatial dynamics. These hydrographic factors affect polyp physiological tolerances and success bay-wide (Cargo and King 1990, Purcell et al. 1999, Shahrestani and Bi 2018. Distribution patterns of C. chesapeakei across sites revealed the greatest concentration of medusae occurred in SLC as opposed to channel sites in UPX and LPX. The inclusion of creek habitat as a source of medusae in our research allowed us to observe an additional strobilation event at the end of the summer from 24 July 2017 to 11 August 2017, when observations of medusae increased from 747 to 4505 medusae in the creek and decreased in the channel from 2048 to 990 medusae. Had we only sampled the channel, the abundance of C. chesapeakei would have appeared to decline through the season, although a new patch of medusae was detectable in the creeks during the same time and formed the second peak in the SLC in August 2017. Although bay nettles are weak swimmers, higher densities of C. chesapeakei in the creek throughout the season suggested that they were able to stay in source habitat and have only a portion of their population transported to the river channel. GAM results showed medusa abundance to be 59 higher in the creek during the peak abundance week of 2016, although volume was calculated as 109 higher in the channel (Table 1). This was not the case during the peak week of 18 July 2017, when a relatively similar number of medusae were found in the creek and channel. The difference in the observed dispersion patterns is not clear, and advective factors such as wind, river discharge, and geomorphology of sites may have resulted in more dispersal out of the source regions (Purcell et al. 2000, Graham et al. 2001, Suchman and Brodeur 2005, Decker et al. 2007, Hamner and Dawson 2009, Kaneshiro-Pineiro and Kimmel 2015, Tay and Hood 2017. However, there were no significant differences in Patuxent River discharge from June to August when comparing 2016 to 2017. Furthermore, summertime wind was not significantly different between 2016 and 2017. Approximately 37% of summer winds blew in from the southwest, with an average speed of 7.8 m À1 s À1 in both years. Similar to the dispersal of mangrove seedlings to shallow coastal habitat via tidal and coastal processes (Duke et al. 1998), medusae need to be dispersed in range of optimal habitat for their planula larvae to settle successfully (Brewer 1991, Lucas 1996, 2001, Sponaugle et al. 2002, Arai 2009). Considering the life span and swimming capabilities of a planula larva (4 d) gives clues to the spatial scale at which physical features such as substrate and flow field, chemical characteristics such as salinity and temperature, and biological processes such as aggregations may be affecting successful settlement within a local habitat (Cargo 1979, Brewer 1984, Brewer and Feingold 1991, Pitt 2000, Lucas 2001, Hamner and Dawson 2009, Holst and Jarms 2010, Lucas et al. 2012, Aglieri et al. 2014, Pitt et al. 2014, Dong et al. 2015. Once planula settles and grows into benthic polyps, physiological cues lead to transverse fission or strobilation (Arai 1997), whereby multiple juvenile jellyfish are released into the local source habitat from a single polyp or strobilae. Further investigation is necessary to understand the mechanisms which control the transport rate and timing of C. chesapeakei movement from source to dispersal habitat.
Previous studies conducted in the Chesapeake Bay mainstem found C. chesapeakei distribution was correlated with patterns of temperature and salinity (Brown et al. 2002, 2013, Decker et al. 2007), which may be due to the coupling between biogeochemical and biophysical processes in the bay (Hagy et al. 2000, Hong andShen 2012). However, we ruled out temperature and salinity effects on the within-season patterns of medusa transport from source to sink, as there were no significant differences in salinity and temperature from the creek to channel habitat within or between sampling events of a given year. For example, on 18 July 2017 we estimated 38 medusae/100 m 3 in the creek and 9 medusae/100 m 3 in the channel. A week later on 24 July 2017, we counted 2 medusae/100 m 3 in the creek and 4 medusae/100 m 3 in the channel, although temperature and salinity were consistent between weeks and sites.
While the mechanisms controlling transport of C. chesapeakei from source to sink were not clear in our study, one potential explanation is density-dependent dispersal. We observed relatively lower abundance and estimated less dispersal in 2016 and higher abundance with more dispersal in 2017. By remaining in the creeks during low-density years, C. chesapeakei lessen the risk of being transported to new, but potentially unfavorable habitats. In a high-density year like 2017, transport to new habitat would be less risky and more favorable as a strategy for population success, as expanding their distribution may potentially mitigate their competition for resources. Density-dependent competition has been postulated as an important mechanism controlling optimal foraging and growth in marine fishes (Fortier and Harris 1989), but it has not been fully explored in jellyfish. Melica et al. (2014) examined density-dependent growth of a population of Aurelia aurita polyps, and Brodeur et al. (2008) reported that jellyfish abundance in the southeastern Bering Sea negatively impacts prey abundance. Further research on the impact of density-dependent factors such as space and food availability on dispersal, growth, and mortality could lead to a better understanding of jellyfish dynamics.

CONCLUSION
A multi-scale approach to ecological studies is often necessary to capture the scale-dependent spatial patterns which likely reflect different controlling factors. For C. chesapeakei, proximity reproduction and turbulence could be important factors for small-scale aggregations and hydrographic conditions and habitat structure could contribute to the difference among different habitats. The distinct difference in patch sizes and dispersal patterns between the July peak period and August peak period in 2017 was likely caused by advective dispersion. Our results suggest that mesoscale spatial patterns of jellyfish might not be caused by salinity, temperature, discharge, or wind events, but instead medusae may have exhibited densitydependent dispersal. The increase in patch size during peak periods and higher estimated dispersal rate when nettle density was high suggest that density-dependent factors such as space and food may play an important role in shaping nettle spatial distribution.
The spatial-temporal patchiness exemplified by our target organism highlights the importance of a research design that captures the life history of a species as it unfolds in space and time, though sampling in shallow habitat is difficult. The multiscale analyses in this study were only feasible due to the collection of high-resolution spatial data, otherwise unattainable with common methods (nets and seines) in marine environments, whereby multiple observations are spatially integrated into a single point. Innovative methods such as sonar imaging proved successful in our research, mitigating the difficulties of sampling in shallow habitats by providing continuous count data along a transect, while providing high spatial resolution to quantify spatial patterns of our target organisms at multiple scales. Understanding the critical role jellyfish species play in a food web exemplifies the importance of interpreting and predicting patterns of scyphozoan medusa dispersal. These methods, in combination with the optimization of mathematical models, will lead to an improved mechanistic understanding of the processes that control species dispersal in marine environments and the potential consequences on the population dynamics of jellyfish species.