Modeling the spatio-temporal distribution of Karenia brevis blooms in the Gulf of Mexico

Study area

This study covered three US states that border the northern GoM, namely the state of Florida, Alabama, and Mississippi (24° to 31°N and 80° to 89°W; see Supplementary Figure 1, Extended data). Supplementary Figure 1 shows study regions in the Gulf of Mexico of K. brevis blooms in the period between 2010-2020 plotted in Google Earth. The identified areas are (1) Southwest Florida, (2) North Central, (3) Central west Florida, (4) Alabama, and (5) Mississippi.

In GoM, the circulation and water inputs are dominated by the Straights of Yucatán, the Mississippi-Atchafalaya system, the Usumacinta River, and the Mobile River.^33,41 The Loop Current, shaped by freshwater inflow from rivers and altered through water density differences and bathymetry, influences the general circulation transport near the Louisiana, Mississippi, and Alabama coastlines.⁴¹ Mesoscale and sub-mesoscale features associated with coastal topography and coastline irregularities can promote primary productivity. These features include upwelling filaments, fronts, shelf and oceanic cyclonic and anticyclonic eddies, and coastal countercurrents.^42,43 West Florida in GoM is subjected to seasonal coastal upwelling, usually intensifying during the early spring to late summer. Along the meridional (southwest) coast, summer northerly winds force offshore transport of near-surface waters, promoting upwelling-favorable conditions.^44,45 During fall, this regime is shifted to southerly winds, favoring downwelling conditions. Along the zonal (south) coast, upwelling conditions are associated with strong westerly winds, but the upwelling intensity is lower. During upwelling relaxation, there is connectivity between the east and west sectors in the Gulf.⁴² Precipitation exhibits a general decreasing trend from the east (Florida) to the west (Texas), while surface evaporation rates generally increase from east to west.⁴⁶

The monthly mean tidal range is low in the winter and high in the late summer, with the highest range exhibited at stations in the central northern regions attributed to monthly changes in seawater density and atmospheric pressure changes.⁴⁷ The salinity in the northern GoM estuaries is influenced by water exchange between the estuarine entrance and the coastal zone and local forcing (tidal advection, river discharge, precipitation) within the estuary.⁴³ River discharge and storm conditions influence the salinity and temperature conditions of offshore waters, making the thermocline deeper and affecting circulation throughout the GoM.⁴⁸

HAB data

HAB cell count was collected from the National Oceanographic and Atmospheric Administration (NOAA) National Centers for Environmental Information’s (NCEI) Harmful Algal BloomS Observing System.^19,49 HABSOS is a system that gathers and distributes data that caters to and provides tracking and observing present and historical harmful algal bloom (HAB), mainly of K. brevis events in the GoM. Data are gathered from the US states such as Florida, Mississippi, Alabama, and other countries being served by NOAA. Florida’s cell counts are monitored and reported by monitoring entities networks; for example, the Florida Fish and Wildlife Research Institute (FWRI), Mote Marine Laboratory (MML), Sarasota County Health Department (SCHD), and Collier County Pollution Control and Prevention Department (CCPCPD). The Mississippi Department of Marine Resources (MDMR) collects monthly phytoplankton samples for routine and responsive (emergency) sampling in 15 stations in the Mississippi Sound. The samples undergo qualitative analysis followed by a quantitative analysis if a toxic bloom is present. Site-specific sampling is conducted by MDCR upon receipt of the bloom report, in addition to trained field technicians’ observations during sample collection upon completion of sample processing. The states submit the data to HABSOS at NCEI, where the data is archived and updated annually.³² In Alabama, HABs are collected and monitored by the Alabama Department of Public Health -Seafood Branch (ADPH) with the Alabama Department of Environmental Management (ADEM) under the regulation and guidance of the United States Food and Drug Administration (FDA) (HABSOS).

In this paper, K. brevis cell count expressed in cells per liter (cells/L) were modeled. Cell count information was collected in the states of Florida, Alabama, and Mississippi, which are the areas in the GoM with the highest presence of HABs in the HABSOS database.

Predictor variables

The environmental variables considered to be potential predictor variables for the habitat modeling were nitrate (NO3), phosphate (PO4), silicate (Si), salinity (Salinity), sea surface height (SSH), mixed layer depth (MLD), sea surface temperature (SST), current velocity (Vel), eddy kinetic energy (Ke), and heading of the current (Heading) (Table 1). Table 1 shows a summary of the environmental variables used in the species distribution models for K. brevis blooms presence (≥100,000 cells/L). All variables were extracted with a 1/4° spatial and daily temporal resolution in the period 2010-2020 from the Copernicus Marine Environment Monitoring Service (CMEMS). Eddy kinetic energy (Ke), current velocity (Vel), and heading (Heading) of the current were derived using the eastward (Uo) and northward (Vo) current vectors. The resolution of the spatial and temporal variables is 0.25° and daily, respectively.

Different environmental covariates were chosen as a proxy of biological, chemical, and physical accumulation in the GoM that can be linked with the presence of HABs. Nitrate, silicate, and phosphate are essential nutrients for the photosynthesis of phytoplankton.^12,39,50,51 These are limiting nutrients, especially in oligotrophic water.^12,39,50,51 Salinity reflects the biological preference of organisms and freshwater input from land; however, in this study, it was used as a proxy for longitude^37,52–54 because of the high correlation between the two (Supplementary Figure 2, Extended data). SSH is associated with mesoscale processes and upwelling.^{13,38,51,54,55} Colder temperatures and MLD are signs of upwelling events that bring nutrients to the surface layer.⁵⁶ Sea surface current flow (current velocity, eddy kinetic energy, and current heading) could also be linked to wind and physical accumulation and transport of HAB cells from the offshore environment to the coastal area.^27,37,52

Environmental variables were extracted from the marine ocean product of the EU Copernicus Marine Environment Monitoring Service (CMEMS). Globally, the Copernicus platform compiles information regarding the physical state, variability, and dynamics of the ocean and marine ecosystems.^55,57 Environmental variables derived from models were selected over in situ observation and satellite products because they offered more spatial and temporal coverage and were less affected by cloud cover. This study used Near Real Time (NRT) and the past reanalysis (RAN) data. Since RAN only covers 2010–2018, NRT data was used for 2019–2020. All environmental variables matched the sampling data (latitude and longitude, observation date).

Modeling HABs presence

A univariate GAM was used to explore the relationship between the presence/absence of HABs (≥100,000 cells/L) in the GoM and different abiotic and biotic variables. The presence/absence of HABS was modeled as follows:

Where g is the link function, which is logit for binomial family, μi is the probability of presence, α is the intercept, f_n are the smooth functions, and x_n are the covariates.⁵⁸

The covariates modeled were the environmental variables (Table 1), spatial variables (Longitude and Latitude), and temporal variables (Month and Year). Month and Year variable were included to account for the seasonal and inter-annual variability in K. brevis blooms, respectively.

A GAM was used in the modeling process to model the presence and absence of K. brevis blooms (≥100,000 cells/L) in GoM to provide a probability of the presence of HABs, irrespective of abundance.⁵⁸ The binomial distribution was chosen because some areas have no presence of HABs, while others have more than 100,000 cells/L. Moreover, areas with abundance equal or higher than the threshold (100,000 cells/L) can have detrimental effects on the biological life in the region and significant socio-economic impacts on people.^59–61

The smooth function degrees of freedom for each explanatory variable were restricted to avoid overfitting. The basis dimension of smoothers was set to k = 6 for the main effect and k = 20 for first-order interaction effects.^62,63 GAMs were fitted using (i) a cyclic cubic regression spline for the variable month, and heading of the current to account for cyclical effects, (ii) a Duchon spline to account for spatial effects (latitude and salinity) considering Euclidean distance, and (iii) thin plate regression splines for all other covariates.⁶⁴ GAMs were fitted using the R package mgcv v1.8.^58,65 A univariate GAMs was done for each predictor variable to give information on the potential functional shape of each predictor variable and their contribution to the deviance explained (percentage of deviance).^57,66

Correlation and collinearity of predictor variables

Two approaches to reduced correlation and collinearity between predictor variables were considered. First, Pearson’s rank correlation was used to examine all predictor covariates to identify highly correlated variables (Supplementary Figure 2A, Extended data). Only one of the correlated variables (Pearson correlation |r| > 0.6) was subjected to the variable selection process at a time and thus produced different model combinations.⁵⁵ Second, the Variance Inflation Factor (VIF) (Supplementary Figure 2B, Extended data) was used to assess multicollinearity among the predictor variables with a threshold value of 3⁶⁷ using the function of from usdm package v 1.1-18.⁶⁸

Model selection

A forward stepwise approach was used for variable selection to select the final model covariates.⁵⁸ Variables were selected if they improved the model’s Akaike’s Information Criterion (AIC) by a value of 2.0 in the selection process.⁶² AIC (AIC = −2l + 2(p + 1)) prefers models that fit the data well and have few parameters.⁶⁹ Moreover, as an additional measure to avoid overfitting in the selection process, only significant variables (p < 0.05) were kept, and non-significant variables were sequentially omitted.⁵⁵ This approach was followed until all variables in the model were significant. The AICs of all best-fitting models were compared, and the lowest AICs were subjected to model validation (Supplementary Table 1, Extended data). Partial effect plots were used to assess the relative contribution of each predictor variable to the HABs presence for the final model with the lowest AICs.

Model validation

The final models were validated using the k-fold cross-validation (k = 5) method to split the training (80%) and testing (20%) data using the k-fold function from the dismo package in R software.⁷⁰ A five-fold cross-validation technique was used and was repeated five times over different random sets to determine the predictive performance of the model.⁷¹ A five-fold cross-validation technique was used to deal with the time-dependence structure of the data.

The model performance of the predicted and observed values were evaluated using a confusion matrix calculated from the PresenceAbsence package v1.1.9 in R.⁷² The area under the receiver-operating curve (AUC), the ‘specificity, the ‘sensitivity’ and the mean True Skill Statistic (TSS) indices were calculated from the confusion matrix.⁷¹ The AUC is threshold-independent and ranges from 0 to 1, an overall measure of accuracy. Values around 0.5 indicate that the prediction is as good as random, and values close to 1 indicate perfect prediction. This index measures the ability of the model to correctly predict where a species is present or absent.^55,73 Specificity is an index depicting the proportion of observed absences predicted correctly, and sensitivity illustrates the proportion of observed presences correctly predicted. The TSS is an alternative measure of accuracy calculated as the sensitivity plus specificity minus 1. TSS values can range from −1 to +1. A TSS of 0 indicates no predictive skill, while a TSS value of +1 indicates perfect agreement and a value below zero (0) indicates no better than random performance.⁷⁴ As opposed to AUC, a threshold-independent approach for model performance, approaches such as sensitivity, specificity, and TSS indices which are threshold-dependent approaches, were also used to determine model performance. To select the threshold, two methods were explored to transform the probabilities into binary predictions (presence/absence)⁷⁵; (i) “sensitivity is equal to specificity” and (ii) “maximization of sensitivity plus specificity.” Ultimately, the “sensitivity is equal to specificity” threshold was used as it gives the most accurate predictions in cases where the dataset has a low prevalence and avoids omissions (false negative) errors.^41,75 The performance scores from the four indices in combination were used to evaluate the predictive performance of the SDMs.

Model predictions

The probability of the presence of HABs in the final model was predicted with a spatial resolution of 0.25 latitude × 0.25 longitude grid cell and monthly temporal resolution between 2010-2020 using the prediction.gam function of the mgcv package v 1.8-36.⁵⁸ The predictions were made using the environmental conditions (significant variables selected in the final model) present in each time step (month between 2010-2020) at a 0.25 × 0.25 spatial resolution. For the model prediction, the monthly predictions between 2010-2020 were averaged to calculate an overall mean prediction map and the standard deviation across the monthly prediction over the 11 years to capture uncertainty. Finally, yearly predictions between 2010-2020 were calculated by averaging the monthly predictions within each year for the probability of the presence of HABs. The ranges for the monthly environmental values for the prediction for the whole study period were checked to ensure that the values of environmental variables were within the environmental envelope used to calibrate the SDMs. About 1.50% of the environmental data used in the model prediction were outside the range of environmental values used to build the species distribution models for HABs.

Model building, evaluation, and predictions were performed using R version 4.0.2.⁷⁶

Spatio-temporal pattern of HABs

This paper studied the presence of HABs (≥100,000 cells/L) in Florida, Alabama, and Mississippi. Because the presence of HABs was mostly in parts of Florida, Florida was divided into three regions¹ Southwest Florida,² North Central Florida, and³ Central West Florida (Supplementary Figure 1, Extended data).

Figure 1 shows the monthly spatial distribution of the presence (≥100,000 cells/L) and absence (<100,000 cells/L) of K. brevis blooms in the Gulf of Mexico from 2010−2020 (shaded 0.25 × 0.25 grids) from the Harmful Algal BloomS Observing System (HABSOS). The white cells mean no sample was collected in the grid cells. Time step 1–12 corresponds to a month in a year such as 1-January, 2-February, 3-March, 4-April, 5-May, 6-June, 7-July, 8-August, 9-September, 10-October, 11-November, 12-December.

Figure 1. Monthly spatial distribution of the presence (≥100,000 cells/L) of HABs.

Most of the presence of HABs was concentrated in Southwest Florida from -81 to -82.5 longitude and 26 to 27.5 latitude. During January, February and December (winter), the highest HAB presence was observed in Southwest Florida (Figure 1). During March-May (spring), the recorded HABs presence was generally the same in all the sampling areas. Spring is the season with the lowest HAB presence with HABs more offshore in May. In June-August (summer), the highest number of HABs presence has been observed in Southwest and North Central Florida, with HAB presence in North Central Florida more offshore than the rest of the season. From September to November (fall), HABs presence recorded was highest in all areas compared to the other season with the highest recorded HAB presence in Southwest Florida. During this time, HABs started to move shoreward in November.

Figure 2 shows yearly spatial distribution of the presence (≥100,000 cells/L) and absence (<100,000 cells/L) of K. brevis blooms in the Gulf of Mexico from 2010−2020 with 0.25 × 0.25 degree grids, from HABSOS. The white cells mean that no sample was collected in that grid cells. Figure 3 shows temporal patterns of the total number of observed presences (≥100,000 cells/L) of K. brevis blooms in the Gulf of Mexico from 2010-2020; Figure 3A shows the total number of presences (≥100,000 cells/L) of HABs sampling by month between 2010-2020 while Figure 3B shows the total number of observed sets with the presence of HABs by month between 2010–2020.

Figure 2. Yearly spatial distribution of the presence (≥100,000 cells/L) of HABs.

Figure 3. Temporal patterns of the total number of (A) sampling and (B) HABs presence (≥100,000 cells/L).

In terms of the year, HAB blooms were irregular and restricted to one state or could occur in two to three states (Figure 2). Blooms in a specific place can occur in consecutive years or every other year. In 2010, HABs occurred in all identified areas, with no area more prominent than others (Figures 2 and 3). In 2011, HABs presence was observed to increase in Southwest Florida offshore the Naples beach. In 2012 and 2013, HABs’ presence expanded from Naples’s beach up to Tampa. In 2014, HABs presence started to increase in North Central Florida in Apalachee Bay, which had the highest HAB presence in the said year. During these years (2010-2014) HABs were observed in Florida and Alabama while no HAB presence nor sampling was done Mississippi. In 2015, there was an increase in HAB presence in Southwest Florida, Central West Florida, Alabama, and Mississippi, with a noticeable decrease in North Central Florida with highest HABs in Tampa and Panama City beach in the Florida Panhandle. In 2016 and 2017, HABs covered only Florida and Alabama, while no presence was recorded in Mississippi. In 2018, HAB presence was highest in all areas. In 2019, only Southwest Florida had a noticeable HAB presence. In 2020, there was a noticeably low number of sampling and presence of HABs in the study area.

The lowest sampling of HABs was found to occur in 2020 (Figures 2A and 3A). The highest HAB monitoring was done from August to November, with increasing monitoring from 2016-2019. The highest level of monitoring regarding the number of samples collected was from September to November, which coincides with the fall season and the months of January and December, corresponding to the winter season. The presence of HABs was found to be highest in January, and October to December (Figures 2B and 3B) and from 2016-2019. Blooms majorly occurred in fall to winter regardless of year or state (Figure 1).

Generalized additive models

After analyzing all possible combinations of covariates (10 candidate models), the final model was chosen based on the lowest AIC values (Supplementary Table 1, Extended data). Because of high correlation and collinearity among several variables, only one variable from the following five groups of covariates was included in the selection process and SDMs at a time (Supplementary Figure 2, Extended data): 1) latitude and longitude; 2) salinity and longitude; 3) nitrate and longitude; 4) nitrate and silicate; 5) nitrate and salinity; 6) nitrate and latitude; and 7) current velocity and eddy kinetic energy. Since latitude and longitude are correlated, we used salinity to replace longitude since salinity is correlated to longitude but not to latitude (Supplementary Figure 2A, Extended data). Therefore, the interaction of latitude and salinity was used to locate grids with a high probability of HABs and answer problems on spatial autocorrelation.

Figure 4 shows the partial effect of the interaction of latitude and salinity in predicting the probability presence (≥100,000 cells/L) of K. brevis blooms in the Gulf of Mexico in the period 2010-2020 in 0.25 × 0.25 grids. White grids mean that there was no sample collected in that area. The interaction between latitude and salinity highlighted three potential areas (1) Southwest Florida, (2) North Central, and (3) Central West Florida, with a higher probability of HABs (Figure 4).

Figure 4. Partial effect of the interaction of latitude and salinity for predicting the probability of HABs.

The results of the final model for K. brevis were summarized and presented (Table 2). The final model for HABs had 23.2% of the total deviance (adjusted r²= 0.123). This model included as predictor variables: the interaction between latitude and salinity, the environmental variables of SSH, silicate, SST, MLD, and heading of the current (Figure 5, Supplementary Animation 1, Extended data), and the temporal variability of month and year (Figure 6).

Figure 5. Partial effects of individual covariate in the final model of HABs in GoM, 2010−2020.

(A) Sea surface height [m]; (B) Silicate [mg/m³]; (C) Sea surface temperature [^oC]; (D) Mix layer depth [m]; (E) Current heading [degrees].

Figure 6. Partial effects of temporal variables (A) Month, (B) Year in the final model of HABs.

Figure 5 shows the partial effects of each individual covariate in the final model for predicting the probability of K. brevis HABs presence (≥100,000 cells/L) in the GoM from 2010−2020. The selected environmental variables were: A) SSH (m), B) Silicate (mmol m^-3), C) SST, D) MLD (m), E) heading of the current (degree). The shaded regions correspond to t standard errors, above and below the estimate of the smooth terms. Short vertical lines located on the x-axis indicate the values at which observations were made. Figure 6 shows temporal variables selected in the final model of K. brevis HABs (≥100,000 cells/L) in the Gulf of Mexico from 2010−2020; The figure shows the partial effects of (A) month and B) year as categorical predictor variables in the final model for predicting the probability of occurrence of K. brevis blooms (≥100,000 cells/L). The dashed lines indicate the 95% confidence interval for the parametric terms, and the black boxes in the x-axis show the distribution of the data.

The Year variable was modeled as a categorical variable that also contributed to explaining the presence of K. brevis blooms. The Year variable was used to account for other variables that were not included in building the candidate models (Figure 6). Lower probabilities of HABs presence were predicted at the beginning of the study period (the year 2010), followed by an increase from 2011-2013 and a gradual decrease in 2014. There was an increase from 2015 to 2016. In 2017 there was an observed decrease in the probability of presence, followed by an increase in 2018. From 2019 there was an observed decrease until 2020. The piecewise construction of each variable for modeling HABs is provided in Supplementary Table 2 (Extended data), and the individual contribution to the percent deviance (Table 3). The time series of HAB probability of presence (Supplementary Animation 2, Extended data) and associated environmental variables used in building the candidate model is shown in Supplementary Figure 3 (Extended data).

Model performance

The predictive performance of the final model for HABs (Table 4) was high, with AUC ≥ 0.8, despite the very low presence of HABs with (≥100,000 cells/L) (6.82%) than absence or HABs with <100,000 cells/L in the data set used.⁷⁷ AUC values of ≥ 0.75 are considered to have good predictive power and can be used for conservation planning.^66,78 The final model had a model sensitivity value of 0.76 which means that the model can identify areas where HABs are present 76% of the time. The model’s specificity was 0.75, which suggests that the model can identify areas where HABs are absent 75% of the time. The TSS score for the HAB species distribution model was 0.51.

Predicted presence of Harmful Algal Blooms

The overall averaged predictions across 2010-2020 identified one important permanent area, (1) Southwest Florida, with higher probabilities of the presence of HABs (Figure 7A). However, the areas in (2) North Central Florida, (3) Central West Florida in the Gulf of Mexico, (4) Alabama, and (5) Mississippi had relatively lower probabilities of HAB presence and high variability based on the standard deviation (Figure 7B). For more detailed daily probabilities of HAB presence with cells ≥100,000 cells/L please see Supplementary Animation 2 (Extended data).

Figure 7. Overall (A) mean prediction of the probability HABs and (B) standard deviation.

The monthly (Figure 8) and yearly (Figure 9) predictions also confirmed the five areas with higher probabilities of the presence of HABs, yet some areas persist over the year, and others are more seasonal. The probability of the presence of HABs is highest during the fall-winter period (November, December-February) in the five identified areas, suggesting a higher probability of presence during this period. In January all three states are prominent, giving Southwest Florida a higher probability of blooms relative to other states. There is also a formation of bands of blooms near the coast and approximately 120 km away from the shore which continue to be present until March. In February, algal growth is more prominent in Southwest Florida while in March, the probability of HABs decreases and there is no area that is more prominent than the other. In April the formation of more offshore blooms starts off the coast of Florida (-81.5 to -86 degrees longitude) and continues to develop and reach its peak in September. In October the offshore bloom was replaced by a more coastal bloom and Southwest Florida became prominent again. In November, the highest probability in all three states was observed. In December, the probability decreases slightly. The seasonal variability of the environmental variable used in building the candidate model is in Supplementary Figure 4 (Extended data).

Figure 8. Overall monthly mean prediction of the probability of HABs (≥100,000 cells/L) in GoM, 2010-2020.

Figure 9. Overall yearly mean prediction of the probability of HABs (≥100,000 cells/L) in GoM, 2010-2020.

The yearly prediction mainly shows only the permanent areas identified in Florida (Figure 9). Because the central West Florida, Mississippi, and Alabama HAB areas occur only in the fall-winter season, these areas were not reflected in the predicted yearly probabilities of the presence of HABs since it only occurs briefly, and yearly maps take the average of monthly predictions in a year. The interannual changes in the probability of the presence of HABs are driven by both the interannual variability and the changes in the final model environmental covariates (Figure 9). The higher probabilities of the presence of HABs were observed to peak in each area in different years and were more variable in Southwest Florida. Higher probabilities of the presence of HABs occurred in Southwest Florida in 2012, 2013, 2015, 2016, 2018, and 2019. The highest probability of HABs was observed in 2016 and 2018. The interannual variability of the environmental variable used in building the candidate model is in Supplementary Figure 5 (Extended data).

Limitations of the study

This study provides insights into the spatio-temporal patterns and environmental preferences of K. brevis blooms. However, the modeling results should be interpreted with caution as they explain a small proportion of the total deviance (23.3%). The low deviance might be explained by the low percentage of data points which has a cell count of ≥100,000 cells/L (presence of HABs) in the HABSOS dataset. Areas below this cell count are considered a no-bloom area. Also, the use of relatively bigger spatial scales of environmental data contributed to the low percentage of deviance, however since this study aimed to observe spatio-temporal patterns in bigger scales and within a 11-year period, we used 0.25 × 0.25 degree resolution data. Despite the low deviance explained, it is typical for studies to model the spatio-temporal distribution to explain a small proportion of the total deviance (10–19 %).^63,66,108 These variables are usually sufficient to explain the target species’ distribution and habitat suitability.¹⁰⁹ This high percent deviance may be due to the quality and quantity of data covering the whole distribution of the data-richer species and the spatio-temporal scales used to conduct the analysis.

The authors acknowledge the biases associated with this data set. Sampling is concentrated if and when there is a bloom, although there are or have been years with more regular monitoring. Moreover, we acknowledge that a 0.25-degree grid is very coarse for identifying small-scale patterns; however, since the study comprised three states, we used a 0.25-degree resolution. Moreover, the said spatial resolution was proven enough to correlate the presence of K. brevis blooms (≥100,000cells/L) between physical/chemical influences on bloom dynamics, especially on the coast of the study areas.

Several ways were suggested to improve and validate the species distribution models for HABs in the GoM. The information and types of data used for modeling species distributions are also growing and rapidly driven by technological improvements in data collection systems.¹¹⁰ For example, satellite applications are increasingly being used to inform and validate species distribution models. There are emerging statistical approaches (e.g., integrated species distribution models) capable of integrating and modeling different types of data sets and getting the most from the data revolution.^110–112 Lastly, future modeling approaches could also explore the other environment and biological preferences of HABs.

In this study, the model residuals had both temporal and spatial autocorrelation due to the seasonal patterns observed in the sampling of HABs and the clustering of reported cell count in spaces. The authors thus suggest accounting for this autocorrelation in future studies.¹¹³