Multidecadal shifts in fish community diversity across a dynamic biogeographic transition zone

A 21‐year fisheries‐independent monitoring dataset was used to explore fish community diversity across a latitudinal gradient to quantify how diversity has changed and relate those changes in diversity to changes in the abiotic environment. Additionally, this study spans a biogeographic transition zone, providing insight into future species assemblages across regions of relatively high species diversity.


| INTRODUC TI ON
Over the past century, increased temperatures have altered sea level, salinity, pH and dissolved oxygen in the ocean (Breitburg et al., 2018;Rhein et al., 2013;Wong et al., 2014). Changes in the abiotic environment have had impacts on habitat suitability, resulting in altered geographic ranges, seasonal activities, migration patterns, abundances and species interactions (McCarty, 2001).
For many marine species, physiology is one of the primary determinants of habitat suitability; many fish are thermal conformers, relying on water temperature to regulate body temperature and related metabolic rates (Clark, Fox, Viner, & Livermore, 2003). As environmental conditions change, fish populations have three general responses: (1) species expand their geographic distribution as environmental conditions become more favourable, (2) species move accordingly while retaining a comparable geographic distribution as favourable conditions shift in location, or (3) species' ranges contract and populations decline, potentially leading to extirpation or extinction as favourable conditions contract in area or disappear (Cheung et al., 2009;Perry, Low, Ellis, & Reynolds, 2005). Together, individual-and species-level responses to biotic and abiotic factors result in changing species distributions and community assemblages.
Current trends of species distribution suggest many marine species are moving from the tropics poleward as minimum temperatures increase at higher latitudes turning formerly inhospi- Warm-water tropical species are expected to move poleward relatively quickly as most tropical species live close to their thermal maximum and must respond more rapidly than species residing in cooler climate regimes (Horta e Costa et al., 2014;Pörtner & Knust, 2007). However, relative to tropical species, many temperate species have broader tolerance limits to varying environmental conditions which could result in a lower rate of movement in response to a changing environment, influencing local extinction rates (Horta e Costa et al., 2014). Varying rates of colonization should result in an increased proportion of warm-water species contributing to diversity within comparatively higher latitude regions undergoing change (Cheung et al., 2009;Estes et al., 2011;Perry et al., 2005). This process of species turnover has been called "tropicalization" (Cheung et al., 2009;Vergés et al., 2016;Wernberg et al., 2013). Evidence of the process of tropicalization has been documented in birds (Thomas & Lennon, 1999), mammals (Hersteinsson & Macdonald, 2016), butterflies (Parmesan et al., 1999), freshwater and marine fishes (Cheung et al., 2009;Perry et al., 2005) and mangroves (Cavanaugh et al., 2014). However, it remains unclear how and whether tropicalization occurs in relatively species-rich systems and, in turn, whether potential novel species assemblages may alter community diversity dynamics.
To address this knowledge gap, here a 21-year dataset is utilized to explore fish community diversity dynamics to gain insight into multidecadal changes in community assemblages and assess factors contributing to those changes. Specific aims of this study are to: (1) examine multiple indices of diversity across a latitudinal gradient to quantify how fish community diversity may be changing in a highly diverse coastal estuary; (2) explore how an area of relatively high biotic change (i.e., a putative biogeographic transition zone) could be utilized to understand and track changes in species assemblages over time; and (3) assess relationships between changes in fish community diversity dynamics and abiotic environment. Together, these aims provide novel insight into how fish communities assemblages are being reorganized in the light of climate-related shifts in species distributions, inform our understanding of how species-rich coastal marine systems respond to climate change, and generate fundamental knowledge regarding community assembly over multidecadal time-scales.

| Study area
Data utilized in this study were generated from an extensive portion of the Indian River Lagoon, FL ( Figure 1). The Indian River Lagoon (IRL) is one of the largest estuaries in the United States, spanning more than 250 km; it is tidally restricted and relatively shallow (average water depth ~1m), being comprised of a mosaic of essential fish habitats including oyster reefs, seagrass beds, mangrove forests and coastal wetlands (Gilmore, 1977). The IRL is composed of three distinct but connected bodies of water that form the broader lagoon system: Mosquito Lagoon, Banana River and the Indian River proper. The biotic community found in the lagoon is comprised of many species found off the Eastern continental shelf of Florida due to the exchange of individuals through five inlets connecting the IRL to the Atlantic Ocean, resulting in the IRL being one of the most diverse estuaries in North America (Gilmore, 1977(Gilmore, , 1995Snelson, 1983). The gradient of environmental factors resulting from the considerable latitudinal extent of the IRL contributes to the relatively high biological diversity. The latitudinal location of the IRL lies at the transition zone between tropical and subtropical or warm temperate species assemblages constructing a putative biogeographic transition zone at approximately 28°N (Gilmore, 1977(Gilmore, , 1995Snelson, 1983).

| Data collection
Samples were collected, and data were generated by the Florida Fish and Wildlife Conservation Commission (FWC) Melbourne-Indian River Field Laboratory as part of the statewide Fisheries-Independent Monitoring (FIM) programme. A description of the random stratified sampling FIM programme can be found at: https ://myfwc.com/resea rch/saltw ater/reef-fish/monit oring/ fim-strat ified-random-sampl ing/). FIM data can be requested through the contact information for the Melbourne laboratory on the following website: https ://myfwc.com/resea rch/about/ conne ct/locat ions/). Fish were sampled monthly with bag seines and a large haul seine using a stratified random-sampling design. Bag seines were 21.3 m long, dragged for 15.5 m and used to collect juvenile and small adult fish (typically < 10 cm) in areas having less than 1.5 m of water. Haul seines were 183 × 3 m, deployed by boat in a rectangular shape along shorelines and on offshore flats and used to collect larger adult fish (Stevens et al., 2016). Fish were identified and enumerated in the field and released, with a subset of samples returned to the laboratory to verify accuracy of identification. At the time of collection, related environmental variables were recorded including temperature, dissolved oxygen, salinity, pH and conductivity using Hydrolab and YSI units.

| Managing data
Raw data for the analyses were provided by FWC. Data and metadata have been archived on the publically available open access Sea Scientific Open Data publication site (https ://www.seanoe.org/).
The data were truncated to the years 1997-2017 and limited to a geographic range from 27.65°N to 28.81°N latitude (approximately 130 km), to provide a more continuous dataset around the putative 28°N biogeographic break. To simplify analyses, mean monthly abundance counts per groupings (explained in more detail in later sections) were calculated for each year during the study period.
Numerical abundance data were 4th-root-transformed to downweight highly abundant species, usually schooling fishes, and give greater weight to mid-range and relatively rare species. Anchoa mitchilli was by far the most numerically abundant species in the dataset (A. mitchilli abundance was roughly five times greater than the second most abundant species, Lucania parva). Analyses were performed with and without A. mitchilli; results of these analyses were similar, so we present results including all species, including A. mitchilli, so as to not exclude any prominent species (Dornelas et al., 2019).
Original data included several gear types including seines, trawls and gill nets. Seines were used for analyses as they had the greatest continuous spatial coverage. The remaining gear types were not included as they were used inconsistently through time and space thereby introducing potential gear biases. Three types of seines were included in analyses, but because larger haul seine catches showed differences in species assemblages when compared to the two small seine catches as well as the physical difference in nets used, these gear types were analysed separately and referred to as "small" versus "large" seines, moving forward. The study area covers three connected bodies of water, cluster analyses were conducted to determine if those bodies of water accounted for differences in species assemblage. Results of the cluster analyses indicated that the species assemblages in the three basins were similar and therefore could be pooled for subsequent analyses, resulting in a more continuous latitudinal gradient (See Appendix A in supporting information).

| Spatial analyses
To analyse the latitudinal gradient, the study area was divided into 29 four-and-a-half-kilometre bins, referred to as "groups" starting from the southernmost point. These groups were used in de- Cluster analysis was performed by utilizing non-metric multidimensional scaling (NMDS) of species assemblage with chronological clustering allowing for visualization of similarity between points in latitudinal order represented by distance on a 2D plane using package "vegan" in R (Oksanen, 2015). Data were normalized, and distances based on Bray-Curtis dissimilarity, which uses abundance F I G U R E 1 Map of sampling area with (a) representing the Indian River Lagoon and the three basins that comprise it. Black lines delineate boundaries of 5-km bins or "groups" and are labelled, and coloured points represent sampling locations within the larger "regions." Inset maps shows (b) the United States with blue star marking sampling location and (c) map of Florida with the sampling extent represented by a blue box. Total sampling area represents 130 km (latitudinally) of the approximately 250-km estuary data among sites for differences. Clusters were constrained to two groups (k = 2) to identify where the latitudinal biogeographic transition zone "best derived break" occurred over the 21-year study period. ANOSIM R statistic was calculated, values approaching 0 represent little difference among groups, while values approaching 1 implies variation among groups (Cook, 2011). Due to extreme cold fronts and repeated algal "superblooms" generating high variability in the time series following 2011, we focused on a subset of the total time series (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) to analyse trends in the biogeographic transition zone. Distances of movement in the biogeographic transition zone were quantified by calculating the straight-line latitudinal distance from one group to another by calculating the change in bins experienced over the sampling period (i.e., 4.45 km/bin X 29 bins).
Indicator species were identified for the species assemblages north and south of the break each year; the R package "indicspecies" determines which species are characteristic of each region considering both abundance and exclusiveness (Cáceres & Legendre, 2009). This package identifies indicator species for the groups of sites being assessed by producing an "indicator value" (IndVal) derived from the exclusivity to the group of sites analysed (A) and number of sites within a group where the species is present (B) (Cáceres & Legendre, 2009). Climate regimes for each species or taxa analysed were assigned using FishBase as it is the most comprehensive source available within the study region (Froese & Pauly, 2019).
A series of alpha-diversity indices were calculated including species richness, Shannon diversity, Pielou's evenness and Simpson diversity. Species richness is the simplest measure of diversity, considering the number of unique species, but is informative with respect to broad presence/absence of species. Shannon diversity builds upon richness by incorporating abundance data, while Simpson diversity, by giving greater weight to common species, deemphasizes the relative value of rare species in favour of species evenness. Pielou's evenness focuses on the proportion of species being examined and not the number of species present. Together, these four indices describe diversity along a spectrum that shifts focus from species richness to evenness, and comparison of these diversity indices can provide insight into how species composition and abundances change spatially and temporally. Using several diversity measures, as opposed to just one allows a more comprehensive understanding of how a species assemblage is changing and allows for more qualitative interpretation.
To understand how diversity has changed temporally and spatially, pairwise beta diversity was calculated using the Jaccard (presence/absence) and Bray-Curtis (abundance) dissimilarities.
Temporal beta diversity was assessed against time zero, defined as mean species assemblages from 1997 to 1999, and compared to three year averages until 2017, resulting in seven time periods.
Spatially, pairwise beta diversity was calculated between pairs of sites moving northward (i.e., 1&2, 2&3 and 3&4). Broadly, pairwise beta diversity describes the dissimilarity of communities between a pair of sites; however, varying components of beta diversity provide further insight into how communities differ. Presence/absence analysis includes overall dissimilarity, turnover quantifies species replacement from one site to another, and nestedness assesses species loss or gain between two sites (Baselga, 2010;Baselga & Orme, 2012). Beta diversity as described by Bray-Curtis dissimilarity can be partitioned to assess the sources of variation in species assemblages and takes species individual abundance into consideration (Baselga, 2013). Specifically, these components describe the patterns in balanced variation in species abundance (hereafter, species balance) and abundance gradient (hereafter, species gradient). With respect to beta-diversity patterns, species balance captures the idea of a particular number of individuals of one species in one location being replaced by a comparable number of individuals of a different species in a different location, while species gradient refers to individuals being lost without substitution when comparing one site to another (sensu Baselga, 2013). Examining these various components of beta diversity provides insight into a dynamic system when combined with measures of the individual site, or alpha diversity.

| Environmental analyses
Multiple linear regressions were run to determine which environmental variables best explained the trends observed in calculated diversity indices. A negative binomial or Gaussian distribution was selected depending on whether the metric was based on count (richness) or continuous (Shannon, Simpson or Pielou's evenness) data.
Conductivity was collinear to salinity and was therefore removed as an explanatory variable. To further examine environmental variables' relationships to species assemblages, an "envfit" test from the R package "vegan" (Oksanen et al., 2018) was performed to visualize which variables accounted for differences in species assemblage. Values of pH ranged from 7 to 9.6 (x = 8.10; SD = 0.25). Salinity (ppt) ranged from 1.5 (ppt) to 44.9 (ppt;

| Spatial analyses-Broad scale
Cluster analyses support the presence of a biogeographic transition zone located within the IRL. Mean species assemblages within the IRL indicate species in the North and Central regions were more similar than the species assemblage of the southern region, the region below the 28°N break (ANOSIM R = .71, NMDS stress = 0.12; Figure 2). One outlier (1997)

| Spatial analyses-Fine scale
Large and small seine data showed changes in the best derived break in species assemblage between the northern and southern region ( Figure 3). Small and large seine catch data exhibited shifts towards higher latitudes over thirteen years of the study (small seine 21.13 ± 11.61 km; large seine 8.51 ± 2.54 km). Large seines had relatively lower latitudinal breaks with all breaks occurring between bins 10 and 13. Small seine catch data produced a stronger northern trend (1.63 km/yr) than large catch data (0.65 km/yr).
Large seine data have higher mean species richness, Shannon diversity and Simpson diversity, and lower Pielou's evenness than small seine data; however, both sized seine datasets follow similar trends temporally and spatially ( Figure 4). Values from large and small seines at the beginning of the sampling period (1997)(1998)(1999)(2000)(2001)(2002) are more similar across regions with the exception of evenness, Overall beta diversity between temporal events was attributed predominantly to the turnover and species balance components of dissimilarity in small seines but were virtually even in large seines.
In small seines species turnover accounted for 73.9% of dissimilarity and 51.0% in large seines ( Figure 5). Species balance accounted for 71.8% of dissimilarity in small seine catch and 52.1% in large seine catch ( Figure 5). In addition, overall pairwise beta diversity between individual 5-km bins was driven by species turnover in presence/absence-based Jaccard's dissimilarity and species balance

| Indicator species
Indicator species were determined for the broad-scale northern Examining climate regimes of the species above and below the break in the biogeographic transition zone, with a focus on tropical species assemblage, the northern and southern regions combined gained 3.12 (±0.92) tropical species over the study period ( Figure 6).
The southern region had higher mean species richness than the north. Additionally, the southern region, defined as the area below the best derived break each year, experienced a faster rate of increase (0.15 species/year) in the number of tropical species inhabiting the region as compared to the northern region (0.03 species/ year), and accounted for the majority of increase in tropical species richness over the entire study area (Figure 6).

| Environmental results
Results of AIC model selection for most diversity indices (species richness, Shannon diversity and Pielou's evenness) for small and large seine catches revealed the combination of all environmental variables (temperature, salinity, pH and dissolved oxygen) to be the model that best described the data (Tables 1 and 2; ΔAIC = 0).
F I G U R E 3 Location of the "best derived break" in species assemblages across a biogeographic transition zone using (a) small seine catch data and (b) large seine catch data. Colours represent the group with warmer red colours representing lower latitude breaks and blue cooler colours representing higher latitude breaks, with group breaks labelled in white. Black dashed line represents linear regression with 95% confidence intervals in grey. These results show a northern shift in the biogeographic transition zone during years of low disturbance and support the idea of novel species assemblages due to varying rates of transition between gear types Simpson diversity from the large seine data follows suit with the previous metrics (Table 2; ΔAIC = 0.0). However, Simpson diversity from the small seine data differs with temperature alone best describing the data (Table 1, ΔAIC = 0.0). In general, the next best performing models were all closely associated with temperature (temperature, month of year, dissolved oxygen), suggesting temperature had a very large influence on fish diversity within the study region (Tables 1 and   2; small gear ΔAIC < 5.0, large gear ΔAIC < 60).
The environmental fit test, which determines which environmental variables best describe the species assemblages present, supported the results of AIC model selection (See Appendix F in supporting information). Temperature followed by dissolved oxygen best described the empirical data of the variables selected; dissolved oxygen was inversely related to temperature. Salinity was the next best determinate, and pH was the least informative variable (See Appendix F in supporting information). Temporally, the southern region in general has higher values of species richness, Shannon diversity and Simpson diversity and lower values of Pielou's evenness. Spatially, higher latitudes have more variance between small and large seine data than lower latitudes

| Biogeographic transition zone
The assumed biogeographic transition zone in the Indian River Lagoon occurs at approximately 28°N (Gilmore, 1977(Gilmore, , 1995Snelson, 1983). The results from analyses pooled across the 21- year sampling period supports the existence of a community transition within 5km (~28.05°N) of this area, while finer spatial scale cluster analysis produced a noisier signal (Figures 2 and 3). The finer scale analyses indicated the location of the biogeographic break fluctuated over several kilometres through time for both gear types, confirming how variable and dynamic the environment and associated biota can be. Despite this variation, both gear types exhibit a northern trend in transition zone location, although with varying rates (Figure 3). The differences in these rates were most likely attributed to gear bias; smaller schooling fish are more likely to be captured in small seines and larger relatively low abundance fish are more likely to be captured in large seines. Attributes of these taxa, such as mobility and residency, could be responsible for the different rates of northern movement and could ultimately lead to novel species assemblages. Similar to this study, other researchers have shown a poleward trend in species distribution shifts (Cheung et al., 2009;Hawkins, Southward, & Genner, 2003;Perry et al., 2005;Vergés et al., 2014), but there are cases with inconsistent or no responses (Chen, Hill, Ohlemuller, Roy, & Thomas, 2011;Kuhn, Lenoir, Piedallu, & Gegout, 2016). Incorporating species-specific differences in response to climate suggests that instead of identifying a distinct biogeographic break in species assemblages, it may be more accurate to identify and consider F I G U R E 5 Temporal beta diversity over twenty-one years compared to time zero. Species assemblages were averaged in three year time increments creating seven temporal bins. Beta diversity was calculated using a) Jaccard's dissimilarity based on presence/absence data compared and b) Bray-Curtis dissimilarity based on abundance data and compare to time zero (1997)(1998)(1999)  The intended purpose of identifying indicator species for the areas above and below the biogeographic break was to identify taxa that could serve as tools or "canaries in the coal mine," to facilitate resource managers identifying when species diversity and biotic communities may be shifting within a given management domain, enable monitoring of those changes and environmental impacts over time and potentially act as a benchmark for mitigation strategy success (Cáceres & Legendre, 2009). However, it is important to note that caution should be expressed when considering indicator species due to dynamic ecosystems, cryptic species and the variation of life stages among species (Zettler et al., 2013). Therefore, rather than holding to a strict definition, we suggest indicator species be considered more broadly as "representative species" that enable a better understanding of differences in species assemblages when comparing among the study regions and as species to explore further when considering potential management applications of the indicator species concept.
In this study, within the northern region, indicator species may have been influenced by the relatively large area of coastal wet- Exploring species assemblages and breaks that occur in the biogeographic transition zone indicate the area has experienced change over the last 21 years. However, there are several factors that must be considered as contributing to a noisier signal. Systems with a high degree of habitat heterogeneity may be comprised of microhabitats allowing fish to survive in areas that would be too extreme without the buffering effect of these potential refugia (Scheffers, Edwards, Diesmos, Williams, & Evans, 2014). In the Indian River Lagoon shallow impoundments located on the grounds of the Kennedy Space Center may function as these types of refugia. There is evidence that climatically, Earth's tropical band is widening and will continue to widen with anthropogenic climate change (Seidel, Fu, Randel, & Reichler, 2008). If this continues, fluctuations and overall change in the fish community assemblage would be expected into the future.

| Indices of diversity
Biodiversity is not homogenous across the Earth and variation exists; efforts to understand the mechanisms driving those differences are increasing, especially through the lens of a changing climate (Alahuhta et al., 2017;Anderson et al., 2011;Koleff, Gaston, & Lennon, 2003;König, Weigelt, & Kreft, 2017;Nekola & White, 1999;Soininen, Lennon, & Hillebrand, 2007;Soininen, McDonald, & Hillebrand, 2007;Tuomisto, 2010aTuomisto, , 2010bViana et al., 2016). Results of species richness, Shannon diversity, Simpson diversity and Pielou's evenness indicate small and large seine data generally exhibit complementary trends, suggesting these types of fish were affected similarly by the factors that drive these metrics. Shannon and Simpson diversity exhibit similar trends implying there were not large changes in common versus rare species, even though the data were transformed to increase the relative weight of rare species. Pielou's evenness values were generally inverse to those of the additional diversity indices; when there was an increase in species, it was less likely those species would be of equal abundance. The combination of these results helps to describe the species sorting processes that influence change in the fish community and provide a baseline for further analysis as sampling continues.
Beta-diversity results in this study show variation among the most important components driving spatial and temporal dissimilarity. Soininen (2018) conducted a meta-analysis assessing the relative importance of beta-diversity components and have found that turnover is usually five times greater than nestedness, in line with previous evidence (Hill, Heino, Thornhill, Ryves, & Wood, 2017;Tisseuil, Leprieur, Grenouillet, Vrac, & Lek, 2012;Viana et al., 2016). Species turnover/balance components are indicative of species sorting due to environmental variables, while nestedness/ gradient components are related to ordered colonization/extinction events (Soininen, 2018). In small seine catch data, the relative importance of turnover was two to three times greater than nestedness, as was the case when comparing the relative importance of species balance and gradient components of abundancebased diversity metrics. However, large seine catch data show turnover/balance and nestedness/gradient components contribute similarly to overall beta diversity. This result implies that larger more mobile fish species are more affected by colonization/extinction events and less affected by environmental sorting than fish species captured in small seines. This further supports the notion of novel species assemblages being formed in the future as these types of fish species respond in different ways. As research continues regarding the mechanisms, putative causes and ultimate consequences of climate change on the loss of biodiversity, these diversity metrics may be important in informing current scientific debate and facilitate a more fundamental understanding of current and future extinction events on a local to global scale (Díaz et al., 2019;Dornelas et al., 2014Dornelas et al., , 2019Urban et al., 2016;Vellend et al., 2013).
Beta diversity decreasing as latitude increases has been identified in plant and animal communities (Qian, 2009;Qian, Badgley, & Fox, 2009;Soininen, Lennon, et al., 2007). In this study the consistent patterns of change in beta diversity found throughout the southern portion of the IRL, with heterogeneity and dissimilarity increasing as one moves further north contradicts the expectation of lower beta diversity moving poleward (See Appendix D in supporting information). This unexpected result suggests additional factors were influencing the observed trends. One of these factors could have been the result of our study region spanning <200km, while these broader trends in diversity play out over larger spatial scales. Additionally, as mentioned above, there are habitat differences across the breadth of the study region that may have a greater influence on beta diversity at the local spatial scale, when compared to broad spatial scale changes in climate.
Another possible factor affecting the trends could be the influence of the subtropical latitude where the broader trends only present themselves in less species-rich temperate systems found at higher latitudes.
When considering temporal beta diversity, generally all components of beta diversity are increasing when compared to time zero ( Figure 5). These results conflict with many reports of biotic homogenization due to changes in habitat and climate, where species assemblages become more similar as opposed to less similar through time (Anibal et al., 2018;Magurran, Dornelas, Moyes, Gotelli, & McGill, 2015;Simberloff et al., 2013). A possible reason for this observed trend is the increased disturbance events including algal "superblooms" that have plagued the sampling area more consistently in recent years. Disturbance may be increasing the dissimilarity of species assemblages and be related to the intermediate disturbance hypothesis where disturbance heavily impacts biodiversity in an area (Connell, 1978). Further exploration of how algal blooms at varying frequencies act as marine disturbances and potentially impact fish community diversity could generate understanding about this trend.

| Environmental drivers
Reports from the International Panel on Climate Change determined the upper 75m of the ocean is warming globally at 0.11°C per decade (Pörtner et al., 2014). Within our study system, broad-scale temperature analysis indicated water temperatures had increased more rapidly (by 0.92°C over the past 21 years), setting the stage for continued shifts in the fish community within the IRL ecosystem as the climate continues to warm. The greater increase in water temperature as compared to the global mean increase is most likely attributed to local environmental conditions; the IRL is a relatively shallow body of water, which increases the influence of changes in air temperature. Dissolved oxygen has been decreasing globally since the middle of the 20th century and is inversely related to temperature, as well as increasing CO 2 concentrations and nutrient inputs (Breitburg et al., 2018). Low-dissolved-oxygen events can contribute to direct and indirect effects on species assemblages and can result in fish kills (Breitburg, Hondorp, Davias, & Diaz, 2008). Many fish kills have been reported recently in the Indian River Lagoon, with exceptionally large kills being documented in March 2016 and August 2018 (Gray, 2016;Cook Pers. Obs.). Further investigations into dissolved oxygen and fish community dynamics could produce useful insights into understanding, predicting and mitigating these events.
Environmental relationships with biota throughout a biogeographic transition zone like the one studied here are less understood in regard to species distributions (Caselle, Kinlan, & Warner, 2010;Selig, Casey, & Bruno, 2010). This study corroborates findings of earlier studies in that excluding the combination of all possible environmental factors together, temperature and its associated variables best describe the diversity of a study region (Clark et al., 2003;Clarke & Gaston, 2006;

ACK N OWLED G EM ENTS
The authors would like to acknowledge Florida Fish and Wildlife Fisheries-Independent Monitoring programme for the extensive data collection. The authors would also like to thank the Marine Ecology and Conservation Lab for their continued help and support throughout this study. Additionally, much appreciation is given to the editors of this paper for giving thoughtful and supportive comments.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw fisheries-independent monitoring data used in this study are publicly available through the Florida Fish and Wildlife Conservation Commission and can be acquired through the contact information for the Melbourne laboratory on the following website: https :// myfwc.com/resea rch/about/ conne ct/locat ions/). Project data and metadata have been archived through Dryad online data repository and are publically available at https ://doi.org/10.5061/dryad. t4b8g thxc. Additionally, the University of Central Florida maintains a publicly available data repository in which published project data will be archived for access and future use (https ://stars.libra ry.ucf.edu/).