The Hauraki Gulf – Tīkapa Moana – Te Moananui-ā-Toi (hereafter the Gulf) is a shallow embayment (maximum depth ∼ 60 m) on the northeast coast of Aotearoa New Zealand (Figure 1). During the austral summer and autumn (warm seasons), the East Auckland Current maintains a southward flow of warm water and onshore winds prevail. In austral winter and spring (cold seasons), offshore winds promote upwelling of cool, nutrient-rich water from the continental shelf, resulting in increased productivity (Chang et al., 2003). Sea surface temperature ranges from 12°C in August/September to 25°C in January/February (Shears and Bowen, 2017).

Observations were opportunistically collected between 2011 and 2020 from the 20 m Dolphin Explorer whale watch vessel operated by Auckland Whale and Dolphin Safari (AWADS). AWADS ran one 4.5 h trip per day during the cold season (June to November) when weather (<20 kts wind) and passenger numbers allowed, and one to two trips per day in the warm season. Average vessel speed was 20 knots while searching for cetaceans and <5 knots during encounters. Between September 2019 and January 2020, we conducted additional surveys from the University of Auckland 15 m vessel RV Hawere on 2 to 3 days per month during favorable weather conditions (Beaufort ≤ 2).

When cetaceans were encountered, at least two observers collected information about location (GPS coordinates), species of cetaceans, seabirds, sharks and/or pinnipeds, predominant behavioral state, and prey (zooplankton or fish) when foraging occurred. Prey type was identified visually or based on the feeding tactic used by Bryde’s whales, which feature a peculiar head-slapping behaviour when feeding on zooplankton (Izadi, 2018). A total of 1,658 multi-species feeding events were recorded: 880 in the warm season and 778 in the cold season.

Data management was conducted in R Studio (Version 1.2.1335) using R software version 4.0.2. We audited the AWADS data to include only observations reporting foraging behavior and created a meta-dataset of all cetacean sightings from 2011 to 2020.

To investigate MSFA composition, we created a species presence-absence matrix including all records from September 2011 to March 2020 featuring two or more predator species where at least one of them was foraging. Rarely sighted cetacean taxa were grouped based on functional similarities (i.e., similar prey and feeding strategies) (Table 1). Due to changes in the AWADS sampling protocol, taxonomic resolution for seabirds was inconsistent, so functionally similar seabird species were also grouped in higher taxonomic units. Petrels and shearwaters, both procellariids, were kept separate as they are functionally different [shearwaters more actively hunting prey, petrels overall more inclined to scavenging, or plankton feeding in the case of storm or diving petrels (Gaskin et al., 2019)].

To identify recurring combinations of predator species or functional groups in MSFAs, we performed hierarchical clustering using Ward’s (1963) linkage method (see also Latasa and Bidigare, 1998; Kramer et al., 2020) implemented in the function hclust (R Core Team, 2013). As a dissimilarity measure for predator associations we used Jaccard distance, commonly employed for presence-absence data (Borcard et al., 2018) implemented in the dist function. To determine the number of clusters we used the NbClust function (package NbClust), which calculates a spectrum of 26 indices for the significance of differences among clusters (Charrad et al., 2014). We set the maximum number of clusters (k) to 10 to obtain a division of meaningful biological interpretation. The most supported choices were a 10-cluster partition (optimizing five indices) and a three-cluster partition optimizing four, namely Marriot (based on optimum stratification; Marriott, 1970), TraceW (trace of the within-cluster scatter matrix, equivalent to the error sum of squares based on within-cluster differences; Milligan and Cooper, 1985; Maulik and Bandyopadhyay, 2002), TrCovW (trace of within-clusters pooled covariance matrix; Milligan and Cooper, 1985), and Ball (average distance of the observations to their cluster centroids; Milligan and Cooper, 1985). For the benefit of interpretability and manageability of the subsequent analysis steps, we adopted the three-cluster partition.

Once the clusters were defined, we calculated the frequencies of the different predator groups in each cluster along with confidence intervals based on the Agresti-Coull algorithm (Agresti and Coull, 1998; Brown et al., 2001). Influential species for between-cluster dissimilarities were identified by SIMPER (similarity percentages) analysis implemented in the package vegan (Oksanen et al., 2019). We then calculated the contribution of each cluster to the overall number of MSFAs by year and season. To ensure a sufficient number of observations per time step, data were pooled by two seasons, warm (December-May) and cold (June-November). Subsequently, we transformed the cluster proportions into time series objects (function ts, package stats) and modeled their behavior over time using autoregressive integrated moving average (ARIMA) modeling (Chatfield, 1995) allowing for seasonal (SARIMA) models, implemented in the function auto.arima (package forecast; Hyndman and Khandakar, 2008). This allowed us to compare all possible models within a range of values for each parameter [namely: lag order of autoregression (p), degree of differencing (d), moving average order (q), lag order of seasonal autoregression (P), degree of seasonal differencing (D) and seasonal moving average order (Q)], and for a maximum overall order (defined as p + q + P + Q), and select the best-fitting model based on AIC, AICc or BIC (i.e., Akaike, corrected Akaike or Bayesian Information Criterion) according to the method chosen (for details on ARIMA model parameters, we invite the reader to refer to Cowpertwait and Metcalfe, 2009). We set the maximum d to 2 to allow for both constant or constantly changing trends, and D to 1 to allow for seasonal differences. All other parameters were left flexible with an overall maximum order of 16 (allowing, for example, for a SARIMA (4,2,4) (4,1,4) [S] model, where S is the seasonal period).

The best model for each cluster identified through the ARIMA approach was then validated by visualizing the autocorrelation of the model residuals in an ACF-plot and the correspondence of sample and model quantiles by a QQ-plot, and by conducting a Ljung-Box test on the residuals (Hassani and Yeganegi, 2019). We also calculated p-values for the autoregressive and moving average terms to test their significance. Finally, the trend component of each time series was extracted with the decompose function and its significance verified through a Cox-Stuart test [function cox.stuart.test, package randtests (Caeiro and Mateus, 2014)].

To analyze the temporal behavior of MSFA diversity, we used a relative “richness index” (the number of functional groups involved in each aggregation; Moore, 2013), as we were interested in changes in diversity over time and the taxonomically inconsistent sampling of seabirds made standard diversity indices unsuitable for our data. After taking the mean of this value for the warm and cold season of each year from 2011 to 2020, we performed ARIMA modeling on the resulting time series.

To test for changes in Bryde’s whale foraging associations, single-species (whales only) and multi-species (whales and associated predators) foraging events were aggregated by season (warm and cold) and year. The proportion of multi-species feeding for each time-step was transformed into a time series object and the best fitting model was selected through ARIMA modeling. We calculated confidence intervals for these proportions based on the Agresti-Coull algorithm (Agresti and Coull, 1998; Brown et al., 2001).

To investigate temporal trends in Bryde’s whale prey preference, we used data from 2011 to 2020 featuring Bryde’s whale foraging (lunging or head slaps). With few sightings of other baleen whales (Table 1), whales not identified to the species level (27 observations) were considered as Bryde’s whales (as per Wiseman, 2008) and included in the analysis.

For records without information on prey type, this was inferred from the presence or absence of indicative species. Obligate fish feeders (dolphins, Australasian gannets, flesh-footed shearwaters, white-fronted terns; Gaskin et al., 2019) engaged in foraging behavior indicate the presence of fish. Zooplankton was inferred from the presence of plankton feeding birds [storm petrels, diving petrels (Pelecanoides urinatrix)] and from the head slapping of Bryde’s whales. When no other species were associated with foraging whales, we categorized this as feeding on zooplankton, as Bryde’s whales feeding on fish are always associated with other piscivorous predators (Rochelle Constantine, personal observation) and they lunge-feed without head slaps (Izadi, 2018). Therefore, for records containing indications of both fish and zooplankton (e.g., Australasian gannets and storm petrels), but not head slapping, we deduced the prey type was fish.

We calculated the percentages of fish and zooplankton feeding for each season and year. For better visualization of temporal patterns, we initially calculated a time series using four seasons and the corresponding ACF. For more robust numbers of observations per time-step, we used two seasons (warm and cold) for the subsequent analysis. An ARIMA model was fitted to the data and model performance evaluated via diagnostic plots (ACF and Q-Q-plot) and the Ljung-Box test. Model coefficients were tested for significance and the trend component was validated by the Cox-Stuart test.

Cetacean species in these MSFAs were mainly common dolphins (observed on 1,435 events; 86.5%) and Bryde’s whales (n = 545; 32.9%), less frequently bottlenose dolphins (Tursiops truncatus), killer whales (Orcinus orca), false killer whales (Pseudorca crassidens), pygmy blue whales (Balaenoptera musculus brevicauda), and sei whales (Balaenoptera borealis) (<2% each; Table 1). The seabirds most frequently recorded were Australasian gannets (82.9%), shearwaters (68.3%), and petrels (19.5%). Terns, little blue penguins (Eudyptula minor), gulls (Laridae), and shags (Phalacrocoracidae) occurred less frequently in MSFAs (<10% each; Table 1). Six records included hammerhead sharks (Sphyrna zygaena) and four included New Zealand fur seals (Arctocephalus forsteri).

Hierarchical clustering assigned 979 records to the first cluster (henceforth Fish1N), 352 to the second (PlankW), and 327 to the third (Fish2W) (Figure 2A). Cluster Fish1N represented almost exclusively fish-feeding events (99.7%; Figure 2B), with a predominance of piscivorous predators such as common dolphins (100%), Australasian gannets (97.8%) and shearwaters (73.3%), and occasionally killer and false killer whales (4.2%), whilst rarely involving Bryde’s whales (1.2%; Figure 2C). Cluster PlankW included both fish (54.5%) and zooplankton feeding events (45.5%) and a variety of predator groups, with a dominance of petrels (84.7%) and Bryde’s whales (58.9%), and occasionally other baleen whales (5.4%) and toothed cetacean (11.4%) species. Of the 206 Bryde’s whale observations in this cluster, 192 (93%) were zooplankton feeding events. Comparing the frequencies of single procellariid species by cluster using the subset of data containing species-specific information on seabirds (2013 to 2020) revealed a dominance of storm and diving petrels (present in 15.1 and 33.8% of the observations in this cluster, respectively). In cluster Fish2W, Bryde’s whales were present in 100% of the observations (92.4% feeding on fish) and common dolphins almost as frequent (99.9%), followed by Australasian gannets (76.5%) and shearwaters (64.2%). Confidence intervals for the estimates of prey type proportion and of predator frequency estimates can be found in Supplementary Tables 1, 2, respectively.

Similarity percentage analysis (SIMPER) identified petrels, common dolphins, Bryde’s whales and Australasian gannets as the most influential species for differences between clusters Fish1N and PlankW, Bryde’s whales and shearwaters for clusters Fish1N and Fish2W and petrels, common dolphins, Australasian gannets and shearwaters for clusters PlankW and Fish2W (Table 2).

The frequencies of all clusters revealed strong temporal variation (Figure 2D). Clusters Fish1N and PlankW became more frequent over the 10-year period, while cluster Fish2W decreased in frequency. Cluster Fish1N predominated throughout, with frequencies between 37 and 73% (except for summer-autumn 2015, 28%). In the time series of the three clusters, cluster Fish2W showed no significant autocorrelation (Supplementary Figure 1). The same models were selected regardless the information criterion used (AIC, AICc or BIC).

The cluster Fish1N time series was identified as a first-order moving average process without seasonality but with a linear trend [ARIMA (0,1,1)]. Diagnostic plots (Supplementary Figure 2) displayed good correspondence of observed and predicted values, with no significant residual autocorrelation. The Cox-Stuart test supported the existence of an increasing trend (Table 3).

The model for the cluster PlankW time series showed a first-order autoregressive process with a non-linear trend [ARIMA (1,2,0)]. The model and the increasing trend were supported by the diagnostic tests (Table 3), despite relatively weak correspondence of sample and model quantiles (Supplementary Figure 2).

For cluster Fish2W, a third-order autoregressive process with a (downward) linear trend [ARIMA (3,1,0)] was detected. Of the three autoregressive terms, only the first was (highly) significant. Diagnostics presented a good fit to the data and the Cox-Stuart test supported the decreasing trend (Table 3; Supplementary Figure 2).

Diversity analysis employed the same MSFA records as used in cluster analysis (n = 1,658). The number of functional groups involved ranged from two to seven, with a relatively constant seasonal mean (range = 2.96 to 3.39 except for summer-autumn 2015 (2.52); Figure 3A). With the low variability of this time series, the auto.arima approach delivered an ARIMA (0,0,0) model, i.e., a constant phenomenon with random oscillations, but no trend or periodicity (Supplementary Figure 3).

Foraging behavior was recorded during 710 (69.9%) out of 1,017 Bryde’s whale observations (which included 27 unidentified baleen whale sightings counted as Bryde’s whales). Other baleen whale species were identified in 95 encounters (Table 1). Bryde’s whales associated with other species in 80.8% (SE = 2.7%; range = 70–96%) of the cases of each seasonal time-step, except for one outlier in winter-spring 2014 (46%).

The time series of the proportion of multi-species foraging events (Figure 3B and Supplementary Table 3) did not show significant autocorrelation (Supplementary Figure 4) and the best-fitting model was again a “null model” [ARIMA (0,0,0)]. Thus, based on our data, there was no evidence for either a trend or seasonality in the tendency of Bryde’s whales to engage in heterospecific foraging events (Table 3) (see discussion). What did change was the cluster they were associated with. Bryde’s whales in the Gulf prevalently fed in association with other predators, and the proportion of observations with Bryde’s whales feeding alone (overall = 19.3%) could be slightly overestimated due to incomplete reporting of birds, seals or sharks present in low numbers and/or for a short time during a foraging event.

To determine Bryde’s whale prey preference we analyzed the same 710 whale feeding records as above, which included 353 (49.7%) fish-feeding and 302 (42.5%) plankton-feeding observations, as well as 55 (7.7%) records with indications of the presence of both fish and zooplankton. Based on the whales’ foraging behavior (Izadi, 2018), we counted these 55 records as fish feeding events increasing the total to 408 (52.2%). The time plot of percent zooplankton feeding over four seasons suggested an upward trend and seasonal oscillations, with increased reliance on plankton during the colder months (Supplementary Figure 6). The latter were supported by the time series showing significant autocorrelation at lag 4 (Supplementary Figure 5A).

For the warm and cold season time series, a SARIMA (0,1,0) (1,0,2) [2] was selected. Coefficients were highly significant (p < 0.0001 for both SMA1 for SMA2). Supported by good diagnostics (Table 3 and Supplementary Figure 5), this model identified the prey preference time series as a seasonal (warm – cold) autoregressive moving-average process with a linear upward trend.

In terms of MSFA clusters, we expected different combinations of predators based on the type of prey; fish and zooplankton attracting piscivorous and planktivorous species, respectively, plus generalist feeders in both cases. Fish-based MSFAs might or might not include the generalist feeders Bryde’s whales, and this was reflected by clusters Fish2W and Fish1N, respectively. Potentially the same holds for aggregations of plankton feeders; however, given the rare presence of other planktivorous cetaceans such as pygmy blue and sei whales in the Gulf (Olson et al., 2015; Barlow et al., 2018), plankton-feeder associations not including Bryde’s whales are composed of seabirds only and thus not included in the data set. We therefore expected to find a single cluster of zooplankton based MSFAs, which is represented by cluster PlankW. As expected, this cluster revealed a very high frequency of petrels, with a dominance of the planktivorous common diving petrels and storm petrels, a high occurrence of Bryde’s whales and occasionally other baleen whales. However, this cluster also included obligate fish feeders like Australasian gannets and common dolphins, and occasionally coastal species (shags and bottlenose dolphins). This overlap may be due to opportunistic feeders such as gulls and terns, or seabirds feeding on zooplankton co-occurring with schooling fish (Zeldis and Willis, 2015), blurring the boundaries between fish- and zooplankton-feeding events. The association of shearwaters and common dolphins is an example of facultative commensalism, with the shearwaters benefitting from the dolphins’ activity but not critically depending on it (Degrati et al., 2014).

Given the competition between common dolphins, Bryde’s whales and gannets for fish, their high co-occurrence in MSFAs (especially in cluster Fish2W) was interesting. Although previous research identified Australasian gannets and common dolphins as predictors of Bryde’s whale presence (Baker and Madon, 2007; Stockin et al., 2009; Purvin, 2014; Izadi, 2018), direct observations of MSFA dynamics show Bryde’s whale feeding lunges often disrupt and terminate the feeding event (Rochelle Constantine, personal observation), as observed elsewhere for other species such as humpback whales (Megaptera novaeangliae) (Haynes et al., 2011).

Our three clusters could easily be generalized to other contexts for comparative studies of MSFA dynamics, provided the predator species are classified in functional groups based on their foraging ecology. Where additional groups are present, integration of additional MSFA types may be necessary. For instance, where two or more baleen whale species are common, different plankton-feeder clusters and corresponding competitive interactions could emerge, and other MSFA types would be expected where foraging behaviors differ strongly from those in the Gulf, such as associations between bottom-feeding gray whales (Eschrichtius robustus) and scavenger birds (Harrison, 1979).

Even though we did not have fine-scale records of the seabird taxa in feeding aggregations available, assigning functional groups enabled robust designation of clusters. Feeding associations comprised a relatively stable number of functional groups. This absence of a trend in diversity was not surprising. Given the trade-off between costs and benefits in heterospecific foraging, it is likely that a maximum number of functional groups exists which can be involved in a feeding aggregation before competitive effects outweigh the benefits of interspecific interactions. MSFAs in the Gulf are not a recent phenomenon and these predators have evolved long-term strategies that exploit the advantages and drawbacks of heterospecific foraging strategies, so this sort of diversity optimum had certainly been reached before our observations started. Nonetheless, in our study the number of groups involved became more constant over time, supporting further stabilization of an equilibrium (Figure 3A).

The majority of Bryde’s whales observed feeding in the Gulf were involved in MSFAs, consistent with observations by Penry et al. (2011) in South Africa. Bryde’s whales accounted for ∼50% of the differences between clusters A and C, implying that whales tended to associate with some species whilst avoiding (or being avoided by) others. The high overall association rate raises the question whether Bryde’s whales engage in MSFAs because they benefit from them, or if it is mainly the other species joining them to increase feeding success or through shared prey preference. Baleen whales rely on multiple sensory systems for foraging at different spatial scales (Torres, 2017). Given the acoustic detection abilities of Bryde’s whales (Constantine et al., 2018) and the overlap in frequency ranges with dolphin and seabird acoustics, the whales may use acoustic cues to detect foraging aggregations. They possibly exploit olfactory cues to detect patches of zooplankton, as found in humpback whales (Bouchard et al., 2019). Surface feeding activity is a useful visual cue for scavenger species and plankton feeding birds like storm and diving petrels, which become commensals of whales by taking small amounts of food without representing relevant competitors (Harrison, 1979). This is supported by the importance of petrels in explaining cluster dissimilarities (Table 2), resulting from the strong association with Bryde’s whales.

By virtue of their generalist feeding habits and year-round presence, Bryde’s whales and potentially other baleen whales may play multiple roles in MSFAs depending on the species they associate with. They are joiners of aggregations initiated by dolphins, and often terminate them by engulfing large amounts of prey (Kahane-Rapport and Goldbogen, 2018). In addition, they likely act as catalysts attracting other species to these events, and as initiators when they are the first to detect prey and attract others through surface-active behaviors. Thus, cetaceans with flexible diet and foraging strategies (i.e., generalist feeders) may be key players in the complex interactions within MSFAs depending on the target prey and on their associates.

Our findings revealed seasonality in Bryde’s whales’ reliance on zooplankton, with a higher frequency of MSFAs feeding on zooplankton in the cold season, when upwelling enhances its availability (Chang et al., 2003; Carroll et al., 2019). Furthermore, 12 years ago fish used to be the predominant prey (Wiseman, 2008), but preference for zooplankton is increasing (see also Carroll et al., 2019). The prey shifts highlight Bryde’s whales’ foraging plasticity both over the long- and short-term including seasonal shifts between fish and zooplankton, within a single population. This may be important for the whales’ tolerance to environmental shifts potentially affecting prey availability (Colbert, 2019; Hauraki Gulf Forum, 2020).

The relative frequencies of the three MSFA types changed over time, with a marked decline of cluster Fish2W characterized by common dolphins feeding on fish in association with Bryde’s whales, whereas the frequencies of clusters Fish1N - mainly seabirds and dolphins feeding on fish, and PlankW - high occurrence of Bryde’s whales feeding on plankton, are increasing (Figure 2D). With the growing reliance of Bryde’s whales on zooplankton, they are progressively abandoning the association with dolphins feeding on fish (cluster Fish2W) in favor of zooplankton (i.e., increasingly contributing to cluster PlankW), leaving dolphins and associated seabirds to feed on fish as shown by the increased frequency of cluster Fish1N. A recent study using DNA metabarcoding to determine zooplankton communities in the Gulf revealed that even though there were seasonal shifts in community structure, preferred whale prey, i.e., copepods, salps and krill-like zooplankton, were present year-round (Carroll et al., 2019). Stocks of pilchards (Sardinops neopilchardus), anchovies (Engraulis australis), and mackerels (Trachurus spp.) are poorly understood (Hauraki Gulf Forum, 2020), and stock assessments are needed to clarify whether the shifts in whale prey preference and in the frequencies of MSFA types are related to prey depletion exacerbating the competition between whales and dolphins and making their foraging association less advantageous for the whales. Whether the whales are foraging on zooplankton as part of a prey availability driven response or as preferred prey remains unknown. Future efforts should aim to identify specific prey via underwater cameras, drone imagery, or the use of eDNA (see Carroll et al., 2019 for prey identification using scat and matched environmental controls), although the high phytoplankton productivity in the Gulf can make visual observations challenging.

Prey consumption by cetaceans in the Gulf is shifting (see also Peters et al., 2020 on common dolphins), which is of interest for future studies of MSFAs. The effective role of seabirds in functional ecosystems within the Gulf requires reliable access to prey, critical for survival and reproduction. Seabirds are ecosystem indicators of nutritional stress (Dunphy et al., 2020), and any impact on prey availability should be considered alongside restoration efforts (Borrelle et al., 2015). Increased competition as a consequence of reduced prey availability is one of the most likely factors contrasting heterospecific foraging associations (Lang and Farine, 2017). With enhanced knowledge of prey dynamics, conservation measures can be adjusted to protect less-flexible foragers, which are more vulnerable to the depletion or redistribution of prey.

Put aside habitat related differences, MSFAs are ubiquitous in marine ecosystems, and the fitness of most participating species depends on these associations to some extent. Focused analyses of the short-term MSFA dynamics will enhance our understanding of the role each species plays in a particular association and reveal behavioral and relational complexities of MSFAs in time and space. To ensure the functioning of these temporary communities and their respective ecosystems, key initiator and catalyst species need to be safeguarded, and maintaining the abundance of prey organisms at their base is equally important. A better understanding of how they change over time will provide important insights on the effects of global and local changes in marine ecosystems.