Loss of suitable ocean habitat and phenological shifts among grouper and snapper spawning aggregations in the Greater Caribbean under climate change

: Phenological shifts have been observed among marine species due to climate change. Modeling changes in fish spawning aggregations (FSAs) under climate change can be useful for adaptive management, because it can allow managers to adjust conservation strategies in the context of specific life history and phenological responses. We modeled effects of climate change on the distribution and phenology of Caribbean FSAs, examining 4 snapper and 4 grouper species. An ecological niche model was used to link FSAs with environmental conditions from remote sensing and project FSA distribution and seasonality under RCP8.5. We found significant differences between groupers and snappers in response to warming. While there was variation among species, groupers experienced slight delays in spawning season, a greater loss of suitable ocean habitat (average loss: 72.75%), and poleward shifts in FSA distribution. Snappers had larger shifts towards earlier phenology, with a smaller loss of suitable ocean habitat (average loss: 24.25%), excluding gray snapper, which gained habitat. Snappers exhibited interspecific variability in latitudinal distribution shifts. Snapper FSAs appeared more resilient to climate change and occupy wider and warmer spawning temperature ranges, while groupers prefer cooler spawning seasons. Consequently, groupers may lose more suitable ocean spawning habitat sooner due to climate change. When comparing species, there were trade-offs among climate change responses in terms of distribution shifts, phenology changes, and declines in habitat suitability. Understanding such trade-offs can help managers prioritize marine protected area locations and determine the optimal timing of seasonal fishing restrictions to protect FSAs vulnerable to fishing pressure in a changing climate.


INTRODUCTION
Climate change has altered environmental conditions on a global scale and prompted changes in latitudinal distribution, depth range, and phenology of fish species in marine environments (Beaugrand et al. 2003, Edwards & Richardson 2004, Poloczanska et al. 2013. In marine systems, climate change has affected the abundance, spatial distribution, and phenology of species from the base of the food web up to higher trophic level predators (Beaugrand et al. 2011). Impacts on tropical fish species are of particu-*Corresponding author: aschr16@ecu.edu

Loss of suitable ocean habitat and phenological
shifts among grouper and snapper spawning aggregations in the Greater Caribbean under climate change lar concern because warming temperatures may result in unprecedented conditions for fishes not experienced anywhere in the world over recent geological conditions (Asch et al. 2018, Reygondeau et al. 2020. Fish spawning aggregations (FSAs) are temporary gatherings of large numbers of conspecific fish that form, often at predictable times and locations, for the sole purpose of reproduction (Erisman et al. 2017). Transient FSAs are a life history phenomenon in which individual fish migrate from within a large catchment area to congregate and spawn in high densities at very specific locations for relatively short periods (i.e. days to weeks; Heyman & Kjerfve 2008, Biggs et al. 2021). In the Greater Caribbean, at least 37 species from 10 fish families are known to form transient FSAs . Fish populations that spawn in large aggregations are highly vulnerable to heavy fishing pressure due to the ability of fisheries to predict and intensively target seasonal spawning locations (Sadovy de Mitcheson & Erisman 2012, Pittman & Heyman 2020. Climate change adds an additional challenge for transient spawning species, which have adapted spawning times to correspond with specific seasonal climatic patterns (Hare et al. 2016). Impacts of climate change are predicted to affect reproductive function of marine fish, with previous studies having identified spawners and embryos as the most temperature-sensitive stages in the life cycle of fishes (Pörtner & Peck 2010, Asch & Erisman 2018, Dahlke et al. 2020. Warmer than optimal temperatures can affect every stage of reproductive development, including spawning, potentially altering the physiology of spawning populations (Alix et al. 2020). Spawning habitat represents a subset of a species' distribution range (see Fig. 1). The timing of migrations and spawning events, or the thermal habitat suitability at the FSA site, may additionally be altered by warming conditions (Asch & Erisman 2018).
Many snapper and groupers (families Lutjanidae and Epinephelidae, respectively) are important species in the Greater Caribbean in terms of both their ecosystem role and as fishery resources (Polovina & Ralston 1987, Arreguin-Sanchez et al. 1996, Amorim et al. 2018. They are managed and harvested as a multi-species complex in the USA, because these families share similarities in life history and ecological characteristics (Farmer et al. 2016, Stevens et al. 2019). Large-bodied species in the snapper−grouper complex are typically long-lived, have late reproductive maturity, and spawn together in transient aggre-gations (Coleman et al. 2000, de la Guardia et al. 2018. These common characteristics make many snapper and grouper species highly vulnerable to heavy fishing pressure and slow to recover if aggregation sites are overfished or extirpated , Sadovy de Mitcheson et al. 2020). There have been substantial fishing impacts to snapper and grouper aggregations documented in the northern Caribbean, including Puerto Rico, the US Virgin Islands, Mexico, Belize, and the Florida Keys, with overall population declines (Claro & Lindeman 2003, Sadovy de Mitcheson & Erisman 2012. The size and structure of aggregations varies among snapper and grouper species, and the characteristics of these aggregations directly influences vulnerability to fishing pressures (Robinson et al. 2015, Biggs et al. 2021. For example, those species or populations that form a few large aggregations during brief periods tend to be more susceptible to rapid fishery declines than those that form many small aggregations over protracted periods (Erisman et al. 2011, Sadovy de Mitcheson & Erisman 2012. While empirical research focusing on climate change influence on fish reproduction is limited, studies on Caribbean fishes have suggested impacts to habitat availability, increased fishing vulnerability, and range shifts correlated with temperature sensitivity and climate change (Fodrie et al. 2010, Maharaj et al. 2018. The goal of this study was to model potential shifts in spawning aggregation distribution as a result of climate change, examining 8 Caribbean reef fish species from the families Lutjanidae and Epinephelidae. Since temperature is a driving factor in species distribution and is directly affected by climate change, we explored how habitat preferences and spawning locations may change with species that spawn across varying temperature ranges. Temperature has been shown to be the predominant effect driving projected distribution changes in Nassau grouper (Asch & Erisman 2018), so we investigated whether this was the case for other reef fishes that form spawning aggregations in the Greater Caribbean. Species with lower thermal habitat preferences during spawning are hypothesized to be more sensitive to climate change and should exhibit larger changes in their distribution, ocean habitat suitability, and spawning phenology. Species with narrow thermal requirements for spawning may need to adjust their seasonality of spawning to a greater extent to stay within the preferred temperature range. With spring spawners, warming temperatures may cause species to spawn earlier in the year. Species spawning in the winter or fall may experience phenological shifts to spawning later in the year as a result of waiting for seasonal temperatures to cool (Pankhurst & Porter 2003).
In the Caribbean, groupers tend to spawn during cooler, winter months, while snappers typically spawn during the warmer spring and summer (Heyman & Kjerfve 2008, SCRFA 2014. We hypothesized that the differences in spawning seasonality and their thermal spawning preferences may result in grouper species experiencing greater changes to their distribution and timing of spawning as a result of climate change, compared to snappers. A multivariate approach was also taken to determine the influences of additional environmental factors on FSAs beyond temperature. Over a centennial scale, we applied an ecological niche modeling approach to compare FSAs of multiple species to assess distributional and phenological shifts, as well as changes in overall ocean habitat suitability, under climate change.

Study area and species
To obtain records of spawning aggregation sites for all of our study species, we used a database developed through collaboration with specialists and con-taining a comprehensive list of known FSA sites , Asch & Erisman 2018. The database included observations and records of FSAs that were verified by direct observation or pubcount (our Fig. 1; see also Heyman et al. 2013. The data aggregate multiple studies done throughout time periods ranging from 1992 to 2011 (Table 1). This database, along with a literature review and consultation from experts on the target species, was used as a reference for this research to identify spawning aggregation sites for each of the study species in the Greater Caribbean within 11.7° to 32.2°N and 61.6° to 92.4°W (Table 1). Geographic descriptions of spawning locations from the database were used to make minor adjustments to coordinates of FSA sites, since data on spawning aggregation locations were frequently rounded for security purposes to prevent exploitation of spawning aggregations by fishers. These adjustments included moving coordinates off land to the nearest stretch of the coastline and inspecting location de scriptions in both the database and the original literature referenced. Latitudinal and longitudinal coordinates were finetuned to match these descriptions from the literature. Any gaps in the data and missing information on spawning months from the database were filled and cross-checked with the scientific literature to confirm  Kobara et al. (2013) and other references included in Table 1. To show that spawning aggregation sites occur across a more limited area of a species range, we also show data in these maps on the range-wide distribution of each species outside of the spawning season. These data are from 2 sources: the Ocean Biodiversity Information System (OBIS; https://obis.org) and Robertson & Van Tassell (2019). The latter reference has had all observations of species validated based on museum collections. OBIS is a more voluminous database, but taxonomic identification in this database is not validated with the same degree of expertise. OBIS also includes observations of larval and juvenile fishes that may have been advected into areas by currents, such as the Gulfstream. Therefore, these 2 datasets are complementary, given these distinctions.  Tables 2 & 3). Lane and gray snapper had fewer recorded FSA observations compared to other species, but they were included in the list of study species so there could be a balanced comparison between the number of species in both families. Nassau grouper FSAs were previously examined by Asch & Erisman (2018) to develop a prototype species distribution model to examine climate change impacts. Data on this species were included herein for comparative purposes to assess more completely how thermal spawning preferences compared across these families of fishes.

Environmental variables
Satellite data on 7 environmental variables were used to assess their effect on the probability of occurrence of spawning aggregations. These variables included sea surface temperature (SST), seasonal SST gradients, geostrophic currents in the east−west and north−south directions, eddy kinetic energy (EKE), chlorophyll a concentration, and vertical velocity of seawater (i.e. upwelling and downwelling). SST was examined because it influences the distribution of many fish species (Pörtner & Peck 2010, Pinsky et al. 2013. Previous research showed SST to be a driving factor for the timing and distribution of spawning aggregations of Nassau grouper (Asch & Erisman 2018). Seasonal SST gradients were examined to determine temperature differences between subsequent months since some species may be triggered to spawn by directional increases or decreases in temperature rather than by an absolute temperature (Wooton & Smith 2014, Asch & Erisman 2018. Chlorophyll a concentration was used as a proxy for biological productivity at spawning aggregation sites, while EKE, vertical velocity, and geostrophic currents are connected to currents that can influence larval fish retention and dispersal, which can in turn influence the suitability of a potential spawning site (Karnauskas et al. 2011, Donahue et al. 2015. The monthly climatology of SST sensed remotely throughout the Caribbean was obtained through the Advanced Very High Resolution Radiometer (AVHRR) Pathfinder, version 5.0 (NODC 2021). Infrared radiometers, such AVHRR, only sense the upper 10−20 μm of the ocean (Nardelli et al. 2005). However, we used both day and nighttime measurements to minimize the effect of thermal stratification that may occur at the ocean surface during the day. Chlorophyll data were taken from the Hermes Glob-Colour dataset version 3.2 (ACRI-ST 2021), which was used to produce a monthly chlorophyll climatology. Vertical velocities were calculated based on data from the QuikScat SeaWinds scatterometer. This data set was downloaded from the National Oceanic and Atmospheric Administration's Environmental Research Division's Data Access Program (NOAA 2021). Geostrophic current anomalies and EKE data were calculated from mean sea-level anomalies (MSLA) from the AVISO satellite altimetry data repository (AVISO 2021). We used a physical oceanographic naming convention where geostrophic flows in the north−south direction are labeled v, while currents flowing in the east−west direction are labeled with u. Chlorophyll concentration and EKE were log 10 -transformed prior to inclusion in modeling due to their log normal distribution. The Reefs at Risk database was used to obtain information on coral reef distribution to constrain projections of spawning occurrence to areas containing reefs (Burke et al. 2011). While the species distribution analysis did not include potential changes in spawning depth, shallow-water coral reef ecosystems are typically limited to a depth of 30 m or less.

Modeling and data analysis
The modeling and data analysis were conducted using Matlab software version R2018a. The Non-Parametric Probabilistic Ecological Niche (NPPEN) model was used to model relationship data on FSA location and timing and information on environmental conditions from satellite data. The NPPEN model was designed to work with presence-only data and is based on a modified version of the Multiple Response Permutation Procedure (Beaugrand et al. 2011). Using a model that handles presence-only data was necessary for this study because there are no confirmed absences of spawning aggregations. A previous study comparing modeling methods of species distribution with presence-only data showed that a technique based on the Mahalanobis distance had one of the best performances when predicting species distribution based on an independent dataset (Tsoar et al. 2007). Since NPPEN is also based on Mahalanobis distance, NPPEN was expected to produce high model skill compared to alternative methods.
NPPEN was used to evaluate ocean habitat suitability for each species and assess what conditions were preferred for reproduction. All possible combinations of environmental variables were used in the NPPEN model initially, and the model with the set of environmental variables that minimized the corrected Akaike information criterion (AIC c ) was selected to use for developing future projections under climate change (Hurvich & Tsai 1989, Asch & Erisman 2018. Running all combinations of the model with the 7 environmental variables produced 128 possible models, including the null model. For species that had an AIC c with the top models separated by values < 2, NPPEN results from each model were used to make future projections, and then the multi-model mean was used in subsequent statistical analysis (Burnham & Anderson 2002). Akaike weights were calculated for each variable for all species to determine the weighted influence of environmental covariates on distribution of FSAs (Burnham & Anderson 2002). Residual deviance explained (D) was used to assess model skill for each species. Deviance was calculated based on the following formula: where df is degrees of freedom, and φ is the dispersion parameter (Quinn & Keough 2002). A φ of 1 was used, since counts of FSAs should follow a Poisson distribution. D was then calculated as the difference between the null model deviance and the deviance of the selected model. Future climate projections were based on the Intergovernmental Panel on Climate Change (IPCC) scenarios using Representative Concentration Pathways (RCPs). We made projections with RCP8.5, which is considered a high-emissions climate scenario with an 8.5 W m −2 change in radiative forcing by 2100 resulting from anthropogenic impacts on climate (IPCC 2013). This RCP was used for analysis since it is the emissions scenario that recent greenhouse gas emissions have most closely tracked (Peters et al. 2013), although there has been a divergence of current emissions from this pathway in the last few years (Hausfather & Peters 2020). An earth system model (ESM) was used to examine how spawning aggregation sites for each species shift under the RCP8.5 climate scenario from 2081 to 2100 and compared against a historical scenario from 1981 to 2000. This analysis considered 20 yr climatologies and conditions were analyzed under both the RCP8.5 and the historical simulation. The principal model used for this study was developed by the NOAA Geophysical Fluid Dynamics Laboratory Earth System Model (GFDL ESM2M) (Dunne et al. 2012(Dunne et al. , 2013. The GFDL ESM2M model was selected due to the moderate equilibrium climate sensitivity compared to the other atmosphere−ocean circulation models included in the 5th Coupled Model Intercomparison Project (CMIP5) ensemble (Cheung et al. 2016). The GFDL ESM2M model had a resolution of 1° at high latitudes, gradually becoming finer scale with a 1/3° latitudinal resolution near the equator (Dunne et al. 2013, Cheung et al. 2016). The physical oceanography component of the GFDL ESM2M model utilized the Modular Ocean Model version 4.1 (Dunne et al. 2012(Dunne et al. , 2013. The marine biogeochemical model used in GFDL ESM2M was Tracers of Ocean Phytoplankton with Allometric Zooplankton, version 2.0, which included a nutrient−phytoplankton−zooplankton− detritus model with 3 phytoplankton functional groups (diazotrophs, small phytoplankton [pico-and nanoplankton], and large phytoplankton [e.g. diatoms]) (Dunne et al. 2013).
To assess projection uncertainty due to the use of different ESMs, we also examined how NPPEN projections differed among 3 different ESMs. These additional analyses were performed only on cubera snapper, which was used as a demonstration species to assess how inter-model differences might affect spawning habitat projections. Cubera snapper was selected as the representative snapper species due to its high number of observations of spawning aggregation sites.
In addition to GFDL ESM2M, we examined the Max Planck Institute MPI-ESM-MR model (Ilyina et al. 2013) and the Institut Pierre Simon Laplace IPSL-CM5A-MR model (Dufresne et al. 2013). These models were selected since, together with GFDL ESM2M, they span the full range of equilibrium climate sensitivity from across the full CMIP5 ensemble. This combination of models has also been used frequently when developing climate projections for other living marine resources and, therefore,lows comparability of our results with those studies (Cheung et al. 2016, Muhling et al. 2017, Asch et al. 2018, Smith et al. 2021). While all 3 models had a 1° latitude/ longitude resolution, each used a separate climate model grid. Environmental data from each model were re-gridded to a common 0.5° grid to allow for greater comparability (Cheung et al. 2016). Average annual NPPEN output from each model was then averaged across grid cells for model comparison. We also assessed the direction of changes in spawning habitat suitability in each grid cell to visualize the degree of agreement among the 3 ESMs. This analysis was completed for both the RCP8.5 and RCP4.5 scenarios to provide an initial assessment of how results might vary depending on the extent of climate mitigation. By the year 2100, RCP4.5 entails a 47% reduction in changes in radiative forcing compared to the high-emissions RCP8.5 scenario.
Model bias correction was used to optimize the comparability between the climate models and satellite observations. This bias correction was based on the monthly mean value of each individual environmental variable throughout the study area. Bias corrections have been shown to be useful in studying climate change impacts on several tropical marine systems (Logan et al. 2014, Matear et al. 2015. This allows for comparisons with observational data and increased model accuracy as bias corrections create a statistical relationship between the original data and modeled values for each covariate and then apply the resulting correction function to the modeled data (McHenry et al. 2019). Monthly mean correction factors were developed for the historical baseline period by comparing climatological values from ESMs and satellite observations, with these correction factors then applied to future periods when making climate change projections.

Statistics
To further analyze the results from future projections of FSA ocean habitat suitability, statistical metrics were calculated for each species and compared across species. Mean latitudinal shift was calculated as the change in weighted average distribution of FSA projections in kilometers per decade for each species. Weights were based on the area of each ESM grid cell. Phenological shifts in month of spawning from historical and future projections indicated the extent to which spawning seasonality could be affected under the impacts of climate change. Central tendency (CT) of the spawning season was calculated as an indicator of phenological change that corresponded to the near center of a species' seasonal distribution. CT can be used to compare skewed seasonal distributions of spawning habitat (Edwards & Richardson 2004). CT was defined based on the following formula: ( 2) where i is the month, and S i is the probability of suitable ocean spawning habitat in month i. These results were expressed as a phenological shift with units of days per decade. An integrated habitat suitability (IHS) score was calculated as a metric of future habitat loss, defined by the loss or gain of habitat based on an array of environmental factors included in the future projections under the RCP8.5 climate scenario. This integral was calculated as a sum across months and locations for probability of spawning habitat to obtain the total ocean habitat suitability across a species range. IHS was defined by the following formula: ( 3) where h(i,j) is the habitat suitability score in model grid cell i and month j, and i max is the total number of grid cells where coral reef habitat occurs.
This was compared between the historical and future scenarios. Change in IHS was measured as a percent change between the 2 scenarios and then converted to an odds ratio to perform a statistical test to assess the hypothesis that there were differences in IHS change across taxonomic groups. For the other 2 metrics measuring climate change impacts (i.e. CT and mean latitudinal shift), an independent, 2-sample t-test was used to examine if differences in sensitivity to climate change occurred between groupers and snappers. The 95% confidence intervals for changes in central tendency, weighted mean latitude, and IHS were calculated using a bootstrap method (Efron & Tibshirani 1998). All sites with NPPEN projections of spawning habitat suitability were resampled with replacement using a uniform distribution. Sample size for each bootstrap iteration was kept the same as in the original analysis. Central tendency, mean latitude, and IHS were calculated for the future and historical periods from this resampled dataset. Spawning habitat suitability probabilities from all months were used for calculating central tendency and IHS, whereas weighted mean latitude was calculated based only on months with documented spawning during the historical period. This process was repeated 500 times with different bootstrap samples. Based on this distribution of samples, we identified the 2.5% and 97.5% quantiles to calculate 95% confidence intervals for each metric.

Comparison between grouper and snapper FSAs
Data presented herein showed that groupers tended to prefer cooler winter months and historically spawned between December and April throughout most of their range, with the exception of Bermuda (Fig. 2). The selected snappers typically spawned later in the spring and summer between March and September, preferring warmer temperatures (Fig. 3, Fig. S1 in the Supplement at www.int-res.com/articles/suppl/ m699p091_supp.pdf). Modeled spawning probability was maximized for groupers between 24 and 28°C (Fig. S1). For seasonal SST gradients, the maxima modeled probability of encountering FSAs fell between −1.5 and 1°C for groupers, with the exception of red hind, for which the model of best fit did not include seasonal SST gradients (Fig. S2). The modeled probability of spawning habitat for the snapper group was maximized at SSTs of 26.5−31°C and between −1 and 2°C for seasonal SST gradients (Figs. S1 & S2). Significant differences between taxonomic groups were seen in CT of phenological shifts and IHS, while the differences between taxonomic groups for mean latitudinal shift were not significant (Table 2). Groupers had a mean spawning phenology shift of 2.8 d decade −1 between historical and future scenarios, while snappers had a shift of −4.7 d decade −1 ( Table 2). Changes in CT were significantly different across taxonomic groups at p < 0.05 (t = 4.6, df = 6). The positive mean phenology change for groupers suggests that these fishes will shift to spawn later in the season, while the negative value for snappers signifies that the spawning season for snappers will move to earlier in the year. The results of the statistical analysis also indicate a difference between snappers and groupers in terms of changes in the proba-bility of suitable ocean spawning habitat (p < 0.05, t = 3.8, df = 6). Groupers appeared more affected by a greater loss of suitable ocean habitat in comparison to snappers. On average, groupers had a poleward shift in distribution of 15.6 km decade −1 compared to a −3.7 km decade −1 equatorward shift among the snappers (Table 3), but the difference between families was not statistically significant (p > 0.05, t = 1.9, df = 6).  Table 4). Those 3 variables had high Akaike weights, between 0.99 and 1.00, suggesting high influence on FSA distribution and a high AIC weight of the se lected model (Table S1). Spawning occurred primarily between January and April each year based on empirical observations, which was consistent with model predictions. Future projections indicated an overall decrease in spawning probability, with a peak during January through April (Figs. 4 & 2c). CT displayed changes of 4.9 d decade −1 between the historical and future time periods, with the positive change in CT indicating a trend towards later spawning (Table 5, Fig. 5). This was the largest shift in phenology among the grouper species. There was a 71.1% projected loss in suitable ocean habitat for spawning based on the IHS score (Table 6, Fig. 6). Mean latitude of FSA sites was projected to shift by 22.5 km decade −1 in the poleward direction (Table 3, Fig. 7).

Black grouper M. bonaci
Four variables had a substantial influence on black grouper FSA projections (Akaike weights between 0.98 and 1.0) and were selected by NPPEN in the best fit model: SST, seasonal SST gradients, and u-and v-geostrophic current anomalies (D = 139.7, df = 99) (Table S1). Both empirical observations and modeled results illustrate that black grouper spawned be tween December through March at all sites, ex cept Bermuda, where spawning occurred between May and November. Future projections indicate that spawning seasonality will continue to peak in De cember through March but will be less seasonally variable (Fig. 2b). Black grouper was projected to shift spawning to later in the season at a rate of 1.8 d decade −1 as a result of its flattened distribution curve for seasonal spawning (Table 5, Fig. 5). Based on the IHS scores from the model, black grouper is projected to lose 69.2% of its potential spawning habitat and have lower probability of suitable ocean spawning habitat overall (

Red hind Epinephelus guttatus
The environmental variables selected in the model of best fit for red hind were SST, v-and u-geostrophic current anomalies, and log 10 EKE (D = 50.0, df = 26). Red hind was the only species in NPPEN to have EKE selected as a variable influencing FSA distribution (Fig. S4). Akaike weights were above 0.99 for SST and v-and ugeostrophic current anomalies, with the Akaike weight for log 10 EKE at 0.79 (Table S1) Fig. 7).

Nassau grouper E. striatus
Primary oceanographic factors affecting Nassau grouper spawning distribution across large spatial scales included SST, seasonal SST gradients, and vgeostrophic current anomalies in the north−south direction (D = 330.3, df = 280). Akaike weights indicated these top 3 variables exerted high influence on FSA distribution (Table S1). Spawning for Nassau grouper typically occurred between De cember and April, apart from Ber muda, where they spawned during June and July. CT of spawning phenology shifted at a rate of 3.1 d decade −1 later in the season between historical and future scenarios ( Table 5, Fig. 5). The 2081−2100 projections showed spawning will primarily occur between January to March (Fig. 2a). The IHS score indicated Nassau grouper may experience an 82% loss in suitable ocean spawning habitat (  (Table 3, Fig. 7). Of the 8 study species examined, Nassau grouper experienced the greatest loss in suitable ocean spawning habitat and the largest mean latitudinal shift (Fig. 4a).

Cubera snapper Lutjanus cyanopterus
The primary oceanographic factors affecting spawning for cubera snapper across a distribution-wide spatial scale were SST, seasonal SST gradients, and geostrophic current anomalies in the north− south direction (D = 34.4, df = 64, Table 7). These 3 environmental variables had Akaike weights above 0.9 (Table S2). Cubera snapper tended to use warm temperatures to spawn, with an average of 29°C (Fig. S1e). With respect to geostrophic velocity, cubera snapper spawn in areas with little-to-no current velocity present (Fig. S5).
Historical models and observations are consistent in indicating that spawning occurred primarily during April through October. Under RCP8.5, future projections indicate spawning will shift earlier in the year to peak between February and June by 2081−2100 (Fig. 3a). CT displayed an earlier shift of spawning season by −5.2 d decade −1 (Table 5, Fig. 5). This reflects approximately a 6 wk change between the historical and future time periods. Cubera snapper experienced the greatest shift in spawning seasonality compared to other species that were modeled. Its IHS score showed an 18.1% decline in spawning habitat, which was the smallest absolute change compared to other species (Table 6, Fig. 6). However, projections still showed a declining probability of spawning habitat occurrence between historical and future periods (Fig. 8). Of all the snapper species, cubera snapper also had the lowest mean latitudinal shift, 1.3 km decade −1 in the poleward direction (

Mutton snapper L. analis
The 4 environmental variables selected in the model for mutton snapper spawning habitat included SST, seasonal SST gradients, and v-and u-geostrophic current anomalies (D = 190.9, df = 88). Akaike weights for all selected variables exerted a high influence on FSA distribution, with values around 1.0 (Table S2). Based on the NPPEN model, historical spawning patterns during the baseline period showed increases in suitable ocean spawning habitat during May through October, although FSA observations did not include October for mutton snapper spawning aggregation formation. Future projections indi-cate a decrease in spawning particularly during June through October, with the maximal amount of spawning projected to occur during December through March by the end of the 21st century (Fig. 3c). This is supported by the CT results, which showed a shift to earlier spawning by −4.8 d decade −1 (Table 5, Fig. 5). Overall seasonality for mutton snapper spawning still encompassed a wide number of months (Fig. 3c) Fig. 7).

Lane snapper L. synagris
The model of best fit for this species had 3 variables, including SST, seasonal SST gradients, and ugeostrophic current anomalies (D = 184.7, df = 57, Table S1). SST, u-, and seasonal SST gradients had high Akaike weights of around 1.0 (Table S2). Comparison of spawning phenology between historical and future periods showed a shift in historical peak spawning from May through September to January through March, as well as a lowered probability of occurrence of suitable ocean spawning habitat (Fig. 3b) Fig. 5). The NPPEN model calculated a 44.2% loss in suitable ocean spawning habitat from the IHS scores between the historical and future periods (Fig. 6). Lane snapper had the greatest equatorward change in mean latitude across all species at a rate of 10.7 km decade −1 (Table 3, Fig. 7).

Gray snapper L. griseus
Gray snapper was the only study species to have a ΔAIC of < 2 across the first 6 models of best fit, so steps were taken to average the NPPEN results of those top models. The first model included SST and vertical velocity (w) as the 2 variables influencing distribution (D = 39.5, df = 14). Results from vertical velocity data indicate gray snapper occurred in areas with slight downwelling (Fig. S6). The remaining models with a ΔAIC < 2 included a combination of 5 variables. The variables included were not only SST and w, but also seasonal SST gradients and u-and v-geostrophic current anomalies (Table 3). Akaike weights for environmental covariates were more variable for gray snapper compared to other species. SST was the primary influence, with an Akaike weight of 1.0, while the remaining selected variables had weights that ranged from approximately 0.4 to 0.6 (Table S2).
Consistent with the historical model, gray snapper was observed to spawn primarily between June through September. The future scenario from the model suggests spawning will become more variable, decreasing during summer and spiking during November−December and April−May (Fig. 3d). This species is projected to develop a bimodal spawning season under the RCP8.5 scenario. Spawning seasonality is projected to shift by −2.9 d decade −1 earlier based on CT scores (  Fig. 6). An equatorward shift in spawning habitat distribution is projected at a rate of 8.5 km decade −1 (Table 3, Fig. 7).

Comparison of models for cubera snapper
Three ESMs were used to assess how choice of model and climate scenario might impact our results. This analysis was performed solely on cubera snap-per to demonstrate the degree of inter-model and scenario variability. Models agreed in the direction of change in spawning probability by the end of the century for a majority of grid cells under both RCP4.5 (71.2% of grid cells) and RCP8.5 (54.6%) (Fig. 9). Additionally, most grid cells with disagreement amongst models in terms of the direction of change were for smaller changes close to zero. For the RCP4.5 scenario, 53.6% of grid cells with disagreements in projected changes had changes less than 0.2 probability. For RCP8.5, 86.3% of disagreements in projected changes occurred in areas with a probability of spawning habitat of 0.2 or less. Overall, the IPSL-CM5A-MR model projected the largest changes in integrated habitat suitability across the study area, while the MPI-ESM-MR model projected the least amount of change in spawning habitat (Fig. S7).

Variables influencing FSAs
For all species studied, spawning habitat was modeled by a combination of SST and a hydrographic variable, including u-or v-geostrophic current anomalies, EKE, and vertical velocity. A previous study analyzing spawning aggregation distribution shifts for Nassau grouper under climate change was replicated for comparison with other species (Asch & Erisman 2018). Our results were consistent with Asch & Erisman (2018), with SST as the most frequent metric selected in the models as influencing species distribution, indicating the importance of temperature in ocean habitat selection and its influence on FSA distribution changes. Our findings were also consistent with previous work showing that seasonal deviations in monthly SST were among the most important influences on the distribution on large-bodied fish species (Mellin et al. 2016). We found that grouper species spawned at similar temperatures (24−28°C). However, Nassau grouper were the most selective in their range with regards to spawning temperature, while red hind was the most variable (Fig. S1). Snappers spawned within similar geographic ranges, but at warmer temperatures (26.5−31°C) compared to groupers (Fig. 2).
Seasonal SST gradients were included in best-fit models for all species, except for red hind, indicating that the change in temperature across seasons was a factor influencing spawning habitat in addition to the effect of absolute temperature. Certain reef fish species spawn during either warming or cooling periods rather than at a specific temperature (Wooton & Smith 2014). Snappers tended to favor warming temperatures with seasonal SST gradients of −1 to 2°C compared to groupers, which most frequently utilized gradients of −1.5 to 1°C, but with high variability between species (Fig. 3).
Four snapper and 3 grouper species had v-geostrophic current anomalies as a metric selected in their best fit model, while 3 grouper and 2 snapper species had u-geostrophic current anomalies selected. For both variables, the probability of suitable ocean spawning habitat was maximized at current anomalies centered around zero, suggesting spawning occurs in areas with slow currents that may result in a greater probability of self-recruitment. Based on genetic evidence and biophysical models, selfrecruitment has been documented to be common among grouper and snapper populations (Paris et al. 2005, Almany et al. 2013, Kough et al. 2016. Alternatively, this may be due to the fact that this study looked at climatologies of geostrophic currents anomalies averaged over a 20 yr period, which resulted in the majority of values for geostrophic current anomalies being close to zero. Due to the use of climatologies, another interpretation of the results is that deviations from typical conditions in currents were not conducive to spawning. EKE was selected as a metric in the model of best fit for both red hind and gray snapper (Fig. S1), while vertical velocity was only detected as an important metric for gray snapper (Fig. S2). EKE and vertical velocity may affect FSA distribution through creation of conditions that influence larval fish feeding, as well as the advection and dispersal of eggs and larvae. In a previous study, drifters released during fish spawning periods were retained in eddies, indicating eddies may act to retain larvae closer to suitable coastal habitat (Heppell et al. 2008). In another study, drifters released at FSA sites moved more quickly offshore than those released at adjacent sites, which suggested that rapid offshore movement lowered predation risk on fish eggs and larvae (Gladstone 2007). However, red hind and gray snapper, the only species to include EKE and vertical velocity in their selected models, had low Akaike weights of 0.29 and 0.14 in their top models, respectively (Tables 4 & 7). This indicates lower confidence for model selection. Compared to other species, gray snapper had a low number of spawning observations, with a sample size of 16, while the sample size for red hind was 30. This may have impacted results. Alternatively, the reduced influence of EKE and vertical velocity on FSA habitat may reflect the spatial scale of our analysis, which focused on changes at the scale of species 107 Fig. 9. Multi-model mean projected changes in spawning habitat suitability for cubera snapper Lutjanus cyanopterus between the future and historical periods for the RCP4.5 and 8.5 scenarios. Colors in the heat map indicate average changes across the GFDL, IPSL, and MPI models. Cross-hatching: locations where all 3 earth system models projected changes in the same direction distribution. Since EKE and vertical velocity exhibit heightened variability at the scale of kilometers to tens of kilometers, finer-scale spatial data may be needed to detect their influence on FSAs.
We initially hypothesized that results may be characterized by greater uncertainty for lane snapper due to the fact that lane snapper is classified as a resident spawner rather than a transient spawner. This could cause its distribution to be less tightly influenced by environmental conditions in the model. Rather than migrating to specific sites with certain environmental conditions, this species stays within its home range for spawning. Resident spawning sites may serve to minimize costs of migration and time exposed to increased predation risk (Donahue et al. 2015). Also, smaller-bodied reef fish may not have the physical and physiological capacity to migrate over long distances, so they may have evolved to adapt to spawn in more variable conditions close to their home range rather than seeking out specific conditions. Despite the expectation that lane snapper FSAs would be less tied to environmental variables due to resident spawning, the percentage of residual deviance explained by the environmentally driven NPPEN model was higher for lane snapper than for any other species. This surprising observation is a subject worthy of future investigation to better understand the underlying mechanisms.
Chlorophyll was not selected in models of best fit for any species, indicating that spawning habitat may not be influenced by biological productivity. This measure alone may not be enough to capture the influence of productivity in the system, or variability in chlorophyll may be too low to impact model results. Biological productivity is potentially important to FSA locations, but it may not be adequately captured by chlorophyll at the spatial scale of this analysis. Different measures of productivity, such as primary production, zooplankton productivity, or export of energy to the benthos, may be more appropriate to include as a variable in future modeling efforts (Stock et al. 2017). Friedland et al. (2012) found that metrics such as particle-export ratio and the ratio of secondary-to-primary production provided greater explanatory power as factors controlling fisheries production than chlorophyll concentration. Additional research could include examining the influence of the chlorophyll maxima at the thermocline on spawning aggregations, in addition to surface chlorophyll.
Overall, spawning habitat suitability for every snapper and grouper species examined herein could be modeled as a function of oceanographic variables. Deviance explained by environmental variables ranged between 39.5 and 330.3 compared to a null model where equal spawning habitat suitability was assigned to all areas with coral reefs. However, the percentage deviance explained by these variables relative to the null model was low to moderate (3.4−20.0%). Also, in Figs. 4 & 8, modeled peak habitat use during the historical period was not always equal to areas with observed aggregations, which implies a level of uncertainty in the projected distribution among certain species. Looking at the magnitude of mismatch of the FSA observations, 1 observation was not predicted well by the models for snapper, while approximately 3 observations were not predicted well for grouper (Figs. 4 & 8).
These factors suggest that additional variables may need to be considered to fully understand spawning dynamics. For example, geomorphology, coral reef condition, lunar phase, and tidal dynamics have been previously shown to affect FSA occurrence among these species , SCRFA 2014. Moreover, since our research focused on large-scale, climatological conditions, it is possible that additional variation in FSA occurrence could be explained if we were able to incorporate in situ data collected synchronously with spawning (Mannocci et al. 2017). However, due to the patchy nature of FSA monitoring programs across Caribbean countries, such environmental data were not uniformly available across the full distribution range of these species. Also, a number of biotic factors likely influence spawning habitat suitability. Population size will determine what percentage of suitable ocean spawning habitat is used (Planque et al. 2007). Since many FSA forming species are currently depleted compared to their carrying capacity, this may have influenced our ability to model the fullest extent of suitable ocean spawning habitat. Interspecific interactions, such as predation risk and availability of suitable prey, may also influence FSA selection (Fernandes et al. 2013, de Araújo et al. 2014. Lastly, some reef fishes exhibit behavioral dynamics where the return to FSA sites might be influenced by the presence of older fish who guide first-time spawners on their migration (NOAA 2013, MacCall et al. 2019.

Phenological shifts in FSA occurrence
When modeling phenology (Figs. 2 & 3), the maximum probability of spawning habitat generally corresponded well with observed spawning phenology (Gokturk 2021). However, low, but non-zero, probabilities of spawning habitat occurred in model output during months when no spawning has historically been observed. This is indicative that there is a non-zero probability of occurrence of the physical oceanographic conditions associated with spawning, but that spawning may be constrained during these months by other physiological, ecological, or evolutionary factors, such as the time needed for oocyte development and vitellogenesis.
There was a clear pattern in phenological shifts between groupers and snappers that fell in line with the original hypotheses that there would be differences in phenological change between taxonomic groups. Groupers shifted phenology slightly later in the season. Groupers spawn during colder months and consequently cannot shift the timing of spawning by very much to track seasonal climate velocity since they are spawning already in the coldest months of the year. This is consistent with previous research suggesting that climate change will cause marine organisms to track the velocity of climate change and shifts in seasonal timing of temperatures (Burrows et al. 2011, Poloczanska et al. 2013. Since marine biodiversity is greater in the tropics, this is of particular conservation concern, as these areas tend to have greater phenological velocities of climate change (Burrows et al. 2011). In contrast to the groupers, snappers were projected to shift spawning earlier in the year, with peak spawning months moving from summer to spring and even late winter, with some variability between species. Gray snapper was the only species to exhibit a bimodal response in phenology, with climate velocities pointing in both directions, meaning shifts in spawning occurred both earlier and later in the year. This difference between how groupers and snappers react to climatic forcing has been observed in other settings. For example, in the Gulf of California, yellow snapper Lutjanus argentiventris and leopard grouper Mycteroperca rosacea exhibited different reactions to El Niño events, with recruitment and landings of grouper augmented during cool La Niña years and snapper recruitment and landings peaking during warm and wet El Niño years (Aburto-Oropeza et al. 2010).
Several snapper species were projected to undergo extremely large changes in reproductive phenology. This includes cubera snapper and lane snapper, which both experienced a 5 mo advancement in their projected peak spawning (Fig. 3). Similarly, projected seasonality of mutton snapper spawning was nearly reversed between historical and future periods (i.e. months with the lowest probability of spawning became months with the highest probability).
However, when viewing these patterns in terms of CT rather peak spawning month, rates of change were more modest (i.e. −4.8 to −5.8 d decade −1 for these 3 species; Table 5), which is just slightly higher than the mean rate of observed change among marine species (Poloczanska et al. 2013). In contrast, the 5 mo change in peak spawning over a centennial scale corresponded to a rate of −15 d decade −1 . These differences reflect the fact that CT is a conservative phenological metric since it accounts for data from all months of the year (Ji et al. 2010). Generally, changes in seasonal extremes (i.e. peak, first, or last occurrence of a species) are subject to quicker phenological change than measures of mean occurrence (Langan et al. 2021).
A key question is whether the rapid shifts in month of peak spawning among snappers is biologically realistic or simply a function of NPPEN assuming that seasonal distributions will track a species' thermal niche. While the gonadosomatic index (i.e. an indicator of preparation for spawning) has been shown to closely track cumulative temperature exposure in many fishes (Ware & Tanasichuk 1989, Lange & Greve 1997, Gillet & Quétin 2006, Neuheimer & MacKenzie 2014, there are likely thermal limits beyond which this relationship would break down. The realism of rapid changes in spawning phenology might also depend upon the source of energy for spawning. Among fish that use recent food intake to provide energy for spawning (i.e. income spawners), there might be a greater capacity and need to rapidly adapt reproductive phenology to environmental changes compared to fishes that utilize stored energy reserves (i.e. capital spawners; Varpe et al. 2009). However, capital spawning may also allow fish with greater flexibility to use stored energy to react to unexpected conditions. Resident spawners, such as lane snapper (Table 1), who do not need to travel long distances to spawning sites, are more likely to be income spawners, whereas species that form transient spawning aggregations are more likely to be capital spawners.
Nonetheless, phenological changes as large as those projected for cubera, lane, and mutton snapper seem feasible, since observed changes of a similar magnitude have been seen among marine species as diverse as copepods, ctenophores, nudibranchs, shrimp, fish, and seabirds (Mackas et al. 1998, Bertram et al. 2001, Philippart et al. 2003, Juanes et al. 2004, Schlüter et al. 2010, Lambert 2013, Langan et al. 2021. Examples in these papers show observed rates of phenological change exceeding a magnitude of 15 d decade −1 (i.e. the modeled rate of change in peak snapper spawning). These examples of extreme phenological change were frequently linked to increases in temperature, with several of case studies coming from climate change hot spots (Lambert 2013, Langan et al. 2021. However, in several cases, additional factors influenced these extreme rates of phenological change (e.g. introduction of invasive species, changes in fisheries management, alteration of predator−prey dynamics; Juanes et al. 2004, Schlüter et al. 2010, Lambert 2013. Among these examples of large phenological changes, several of them negatively affected the survival, abundance, or recruitment of protected or commercially valuable species that interacted with the organism undergoing the sizable phenological change (Bertram et al. 2001, Philippart et al. 2003, Mackas et al. 2007).

Latitudinal shifts in FSA sites
Compared to phenological changes where there were distinct differences among fish families, the results of this study are less clear-cut when looking at projected latitudinal shifts in groupers and snappers. Results for latitudinal shifts differed among species and ranged from −11 to 40 km decade −1 (Table 3). Groupers tended to have larger poleward shifts in distribution, but this also varied from species to species. Similar to the pattern seen with phenological shifts, Nassau and yellowfin grouper had the greatest latitudinal distribution shifts within the groupers. Equatorward shifts were observed with lane and gray snapper, whereas red hind was projected to have an extremely small shift in distribution (0.01 km decade −1 ) ( Table 3). Models from 5 of the 8 study species that all had high confidence based on Akaike weights showed a poleward shift in distribution. This is consistent with the global trend towards poleward distribution shifts observed for many fish species (Cheung et al. 2009, Morley et al. 2018. Local extinctions and invasions can be caused by shifts in latitudinal ranges (Parmesan 2006, Cheung et al. 2009). These shifts in distribution have generally extended in the poleward direction with increasing temperatures, resulting in greater potential for climate-induced invasions and higher invasion intensity at higher latitudes (Hiddink & Hofstede 2008, Cheung et al. 2009, Fodrie et al. 2010. Conditions that exceed the temperature range of fish species may make reproduction at preferred FSA sites and seasons no longer possible, forcing them to adapt or shift spawning to cooler regions and seasons to avoid extinction (Dahlke et al. 2020).
While this was the first study to examine climate change impacts on multi-species spawning aggregation locations of Greater Caribbean groupers and snappers, several previous studies have investigated distribution shifts among these species outside of their spawning season. Morley et al. (2018) showed less agreement in distribution shift direction of species with low certainty in some models, which was the case for gray snapper, lane snapper, and red hind. The southward shift for gray snapper identified in this study contrasts with previous research indicating a poleward distribution shift (Hare et al. 2012). However, the magnitude of distribution shifts for FSAs in general and gray snapper in particular are dependent on the climate change scenarios examined in each study. Also, the smaller sample size for gray snapper FSA sites and low certainty in model selection may have influenced results. Additionally, this study limited spawning observations to sites with presence of current-day coral reefs and only examined the spawning life history stage, which differed from Hare et al. (2012). Since our study modeled gray snapper and other species as obligate reef spawners and did not make projections of future changes in coral reef distribution, this likely limited the extent of poleward distribution shifts. However, gray snapper can use coral reefs facultatively during spawning and may alternatively use other structured, hardbottom substrate (SEDAR 2018). Future work should explore how this affects projected FSA distribution by parameterizing NPPEN to make projections of spawning in areas with coral reefs, other natural hard substrates (e.g. rock reefs), and artificial substrates, such as oil rigs and pipelines, that can serve as artificial reefs (Paxton et al. 2020).

Loss in suitable FSA ocean habitat
Our integrated habitat suitability metric explained gains and losses in marine spawning habitat independent of shifts in fish distribution and changes in phenology. Overall loss of suitable FSA habitat between historical and future periods under the RCP8.5 scenario ranged between 68 and 82% for groupers (Table 6). This contrasts with the smaller 18−44% loss of habitat among lane snapper, mutton snapper, and cubera snapper, producing significant differences between these 2 groups of fishes (Table 2). Gr ay snapper differed from all other study species and showed a gain of suitable ocean spawning habitat of nearly 105% by 2100 (Table 6). Previous work on gray snapper has also found that this species is projected to ex-pand its range under the impacts of climate change as temperatures increase (Hare et al. 2012, Morley et al. 2018. For example, a study conducted on gray snapper in the Gulf of Mexico found a 71% increase in thermal habitat (Morley et al. 2018). Range expansions may increase local biodiversity on short-term time scales as poleward-retreating species are outpaced by poleward-advancing species. However, in the long term, these climate range expansions may act similarly to nonnative species invasions, modifying local dynamics related to competition, predation, herbivory, and parasitism (Fodrie et al. 2010). Populations shifting to favorable environmental conditions may still experience novel selection pressures from altered biotic interactions and unprecedented combinations of photoperiod cues and climatic effects (Moran & Alexander 2014).
Some research has suggested that rather than range shifts, it may be possible for fishes to respond to climate change by adjusting thermal range through individual acclimatization or evolutionary adaptation across generations (Angilletta 2009, Dahlke et al. 2020). This could be applicable to the study species in their ability to adapt to changing temperatures and spawn in warmer conditions rather than shifting their spawning habitat and range. The changes in suitable ocean habitat for FSAs may also be tied to other ecological needs for spawning habitat in addition to temperature requirements. Spawning sites may provide necessary substrates for egg deposition or hydrographic features that assist with egg and larval dispersal, including current speeds and flow direction (Heyman & Kjerfve 2008, Wooton & Smith 2014, Dahlke et al. 2020. As a result, even if a range shift would provide more suitable temperature conditions for species forming FSAs, additional requirements needed for reproduction may be missing, including factors related to hydrographic conditions, reef geomorphology, and biotic conditions like predator absence and food availability. Habitat fragmentation as a result of range expansion could result in decreased larval connectivity (Moran & Alexander 2014). Studies using coupled physical−biological models suggest larval transport could influence the range of marine species irrespective of local environmental conditions like temperature (Gaylord & Gaines 2000, Thompson et al. 2018. Changes and shifts in spawning habitat have further ecological implications, since predation, competition, and prey availability can change in response to co-occurrence of temperate, subtropical, and tropical species (Sax et al. 2007, Fodrie et al. 2010). These ecological interactions may influence whether it is possible to establish new FSAs in locations with suitable climatic conditions.

Future research
Further research on transient FSAs is needed, as there may be sites that are unmanaged or undiscovered currently ). Knowledge gaps on FSAs limit the effectiveness of management strategies meant to protect FSA habitat and productive fisheries, consequently limiting the performance or design of marine protected areas (MPAs) (Sale et al. 2005, Crowder & Norse 2008. Fish populations can be sustained effectively through no-take reserves if spawning occurs within the boundaries or there is connectivity between FSA sites and marine reserves (Sale et al. 2010).
As establishment of MPAs and implementation of fishing restrictions can be depth-dependent, change in depth of spawning should be considered for future research. Similar research conducted has projected species to become restricted to deeper habitats based on historic observations (Pinsky et al. 2013, Kleisner et al. 2016, Morley et al. 2018). This could be an important metric to consider for species that may adjust depth in response to warming temperatures rather than adjusting their phenology or latitudinal distribution as was the focus of this study.

Conclusions
Researchers should continue to identify and monitor transient FSA sites following standardized protocols and sharing results via a cooperative FSA monitoring network (Heyman et al. 2019, Pittman & Heyman 2020. Future work exploring the life history characteristics of different species may provide key insights into responses under the impacts of climate change, and consequently managers may need to adjust current management strategies. Different life stages of fish species are not equally sensitive to temperature and the effects of climate change. These differences make it critical to develop modeling approaches that consider the interactions between multiple life history stages under the same framework.
Overall, there was evidence from this study supporting that climate change will impact the phenology, distribution, and ocean habitat suitability of reef-fish transient FSAs. Among these effects, strong differences were seen between the grouper and snapper families, with varying life history characteristics and thermal preferences for spawning. For each species, SST coupled with a hydrographic variable was found to influence distribution and proba-bility of spawning habitat. With cooler preferences for spawning habitat, groupers were found to be more impacted by climate change under the RCP8.5 scenario compared to snapper species, which was consistent with our hypotheses. While the directionality and average latitudinal shift varied across species with no significant differences between families, there were differences between families in terms of changes in suitable ocean habitat and phenological shifts. Snappers are expected to shift spawning earlier in the season, while groupers' spawning season was projected to shift slightly later in the year and have a greater loss of suitable ocean spawning habitat. Differences in our results between shifts in mean latitudinal shift, habitat suitability, and phenology highlight the importance of looking at multiple metrics when studying fisheries under climate change. While there are many studies exploring distributional changes, other metrics, such as phenological shifts or overall habitat availability, may provide additional insight into climate change responses. This is especially important since our results suggest that there may be trade-offs between whether species respond to climate change by altering their phenology, mean latitudinal distribution, or by experiencing declines in habitat suitability. The information obtained from this study could be a useful tool in future management of both marine reserves and other types of protected areas, as the timing and locations of spawning events could be altered as a result of climate change. This has implications for harvest and fishing restrictions, as well as seasonal sales bans and site closures during spawning events.