To assess the relative impact of oil on seabirds in the nGoM, we created an index of the vulnerability of seabirds to oiling (i.e., the vulnerability of seabirds to oiling index, VSOI). We developed the VSOI by adapting recently described frameworks and index-based approaches that assessed the potential impacts of offshore energy development on marine fauna and then tailored it to the current state of knowledge of seabirds in the nGoM. Our approach was informed by similar efforts (e.g., Kelsey et al., 2018) and proposed frameworks (Polidoro et al., 2021) and the terms used in our approach are most closely related to those described in Nevalainen et al. (2019). The framework applied by Nevalainen et al. (2019) best reflects the current state of knowledge for seabirds in the nGoM. In brief, seven variables were used to describe the vulnerability of seabirds to oiling. Variables are conceptually grouped into two sub-indices; exposure potential (n = 4 variables) and sensitivity (n = 3 variables). Exposure potential characterizes the probability of an individual co-occurring or potentially interacting with oil, and was described by seasonal occurrence, overlap with O&G platforms, flocking behavior, and primary foraging technique. Sensitivity characterizes the relative impact of mortality on the population and was described by age at first breeding, duration of the incubation and fledging period, and residency status within the nGoM ( Figure 1 ). Section 2.3 describes each variable in detail. We also describe the uncertainty of the information used to calculate the VSOI to identify knowledge gaps (Section 2.4). Lastly, using spatial layers produced to characterize species-specific habitat, we also describe the spatial extent of the habitat of all seabird species (species habitat footprint), the sum of seabird vulnerability (cumulative seabird vulnerability), and the overlap of O&G platforms with cumulative seabird vulnerability (Section 2.5).

We gathered information from multiple sources (described below) and used the information collected to assign a qualitative score to each variable.

The spatial extent of analysis is constrained to the survey footprint of GoMMAPPS ( Figure 2 ). Therefore, the rank of seabirds relative to each other is constrained to this same spatial footprint. Additional habitat for these species may occur within the Gulf of Mexico beyond the boundaries of the study area, including coastal areas that were excluded due to cruise logistics, and the southern Gulf of Mexico, which was not surveyed as part of GoMMAPPS (i.e., south of the U.S. Exclusive Economic Zone, which formed the southern extent of the survey footprint).

To describe the vulnerability of seabirds to oil through the VSOI, we characterized seven variables for each of the 24 focal species ( Table 1 ). These variables were informed through a literature review, GoMMAPPS vessel survey data, and data from the BOEM data center (above). Seven variables are conceptually grouped into sub-indices ( Figure 1 ). For each variable, the relative potential (exposure) for or intensity (sensitivity) of an interaction with oil is scored as 0 (none; if applicable), 1 (low), 2 (medium), or 3 (high). For variables where only two states were possible, scores were either low (1.67) or high (2.33). Scores ( Table 1 ) for each variable were separated either by (1) intuitive intervals (i.e., seasonal occurrence (number of seasons, defined below), overlap with O&G platforms (even intervals), flocking behavior (log scale), residency (binomial), foraging technique (greater or lesser exposure to oil on the water), or (2) observed gaps in the spectrum of values (i.e., age at first breeding, days from incubation to fledge). The sum of the scores for each variable was then used to define the vulnerability of each species to oiling. The VSOI scores characterize each focal species’ relative vulnerability to oil but do not provide or predict the absolute vulnerability to oil.

To acknowledge the imperfect nature of the data used to score each variable and the natural variation associated with each variable, we assessed the extent to which uncertainty can impact our estimate of vulnerability. Here, we consider uncertainty to encompass a lack of species-specific information, low or poor rigor, and natural variation. Based on the number and quality of data sources available and the level of natural variation, we assigned quantitative values for the characterization of the uncertainty of each variable for each species as either low (10%, 0.10), medium (20%, 0.20), or high (50%, 0.50), e.g., Kelsey et al. (2018). For example, high uncertainty (0.5) would occur if no data were available for the age at first breeding for a given species and values were obtained from a related species. In contrast, low uncertainty (0.1) would occur when data were directly available from the study area and where such data could be assumed to be relatively consistent across years. For example, flocking behavior was measured directly during vessel surveys and is presumed to be a behavioral characteristic that is relatively consistent across time. Medium uncertainty (0.2) could include contemporary observations but of a timescale unlikely to capture the full extent of natural variability. For example, seasonal occurrence in the study area was observed, but observations were uneven across seasons ( Figures S1, S2 ).

Calculating the range in the score for a particular variable, including uncertainty, required knowing the: 1) score of the variable, 2) uncertainty for that variable, and 3) difference between the maximum and minimum potential scores for that variable. The range was calculated as the score of the variable +/- the uncertainty multiplied by the range of potential scores. For example, for age at first breeding, a given species has a score of 2, high uncertainty: 0.5, and potential scores of 1 to 3, with the difference being 3-1 = 2. The range of the score for age at first breeding would be 2 +/- (0.5*2); therefore, 1 (lower limit) to 3 (upper limit). The resulting lower and upper limits describe the uncertainty in each variable and, when combined across variables, describe the uncertainty in the vulnerability score for each species. We constrained the lower and upper limits to the extent of the variable scoring range, either 0-3 or 1-3 ( Table 1 ). To identify knowledge gaps, we compared the average uncertainty of each variable where greater uncertainty indicates a greater knowledge gap.

To define the spatial extent of all seabird habitats within the study area (hereafter seabird habitat footprint), we overlayed highly suitable habitat (Section 2.3.1. Exposure potential) of all species combined and used the entire area that was highly suitable habitat for >1 seabird species. This synthesizes species-specific habitats to provide an assemblage-wide seabird habitat footprint. To identify the areas of the highest vulnerability for all seabirds, we assigned the vulnerability score of each species across its high suitability habitat, then summed the areas of overlap between species. This defined the cumulative seabird vulnerability in a spatially-explicit context. The cumulative seabird vulnerability within the O&G footprint was compared to the cumulative seabird vulnerability outside of the O&G footprint and described as the cumulative seabird vulnerability/km². This is a proxy for the intensity of seabird vulnerability in each area. Thus, the cumulative seabird vulnerability was used to characterize the distribution and intensity of the vulnerability of the seabird assemblage in the nGoM to oiling.

To review our characterization of high suitability habitat overlapping O&G platforms, seabird habitat footprint, and cumulative seabird vulnerability: 1) the habitat overlapping with O&G platforms is a variable that defines the ‘spatial extent of potential exposure to oil’ of an individual seabird species as a percentage of highly suitable habitat, 2) the seabird habitat footprint utilizes highly suitable habitat and characterizes the area within the ‘spatial extent of all seabird habitats’, and 3) the cumulative seabird vulnerability synthesizes highly suitable habitat and vulnerability scores to describe the ‘distribution and intensity of the vulnerability of the seabird assemblage in the nGoM to oiling’. Unless stated otherwise, all other data analysis was performed in R version 4.1.1 (R Development Core Team, 2021).

Scores for vulnerability to oiling for seabirds in the nGoM ranged from 10.33 - 19.33 (median = 14.67, mean = 14.47 ± 0.44 SE), with a potential range of 5.67 - 20.33 ( Table 2 ; Figure 3 ). Species in the lower 20% of VSOI scores were: Parasitic jaeger (Stercorarius parasiticus; 10.33), Bonaparte’s gull (Chroicocephalus philadelphia) and Pomarine jaeger (Stercorarius pomarinus; 11.33), and Brown noddy (Anous stolidus) and Common tern (Sterna hirundo; 12.33). All but Brown noddy, which breeds in the Dry Tortugas and on the Campeche Bank in the southern Gulf of Mexico, breed in the continental interior of North America or the high Arctic. The five species with these lowest vulnerability scores include only two families: Stercorarius and Laridae. The relatively low scores for these species were due to low values across four variables ( Table 3 ). The incubation through fledge period is relatively short in these species compared to other species, with four of the five species spending less than 20% of the year with a chick ( Table 2 ). The primary foraging techniques (i.e., kleptoparasitism and surface feeding) resulted in less exposure than foraging techniques of other species (e.g., plunge-diving). The spatial overlap between the highly suitable habitat and O&G platform footprint was relatively low compared to other species. Lastly, the age at first breeding for these species was relatively low (≤ 4 years) compared to the average age at first breeding of other species ( Table 3 ).

Species in the upper 20% of VSOI scores were Northern gannet (Morus bassanus; 19.33), Audubon’s shearwater (Puffinus lherminieri; 18.33), then Brown booby (Sula leucogaster), Great shearwater (Ardenna gravis), and Herring gull (Larus argentatus; 16.33; Table 2 ; Figure 3 ). Breeding ranges for these species include maritime Canada, the interior of North America, the Caribbean, and the South Atlantic. These species include three families: Sulidae (Northern gannet, Brown booby), Procillariidae (Audubon’s shearwater, Great shearwater), and Laridae (Herring gull). The relatively high scores for these species were due primarily to high values across three variables ( Table 3 ). The primary foraging techniques (i.e., plunge- and pursuit-diving) resulted in higher exposure than foraging techniques of other species (e.g., dipping). The incubation through fledge duration was relatively long in these species, with four of five being > 30% of the year, compared to the majority of the other species (i.e., 20-30% of the year). Lastly, although there was variation in the age at first breeding, the age at first breeding for these species tended to be higher than the remaining species.

Three dominant rank-forming variables were shared in separating the species with the least and greatest vulnerability ( Table 3 ). Two were related to sensitivity; breeding age and days from incubation to fledge, and one was related to exposure; primary foraging technique. The percent of habitat overlapping O&G platforms was important when separating species with the least vulnerability from the other species but did not have a notable impact on separating the most vulnerable species from the other species.

The average uncertainty among species appeared greatest for sensitivity variables compared to exposure variables. The uncertainty surrounding age at first breeding (0.38), number of days incubation through fledge (0.37), and residency (0.24) were relatively high. In contrast, the uncertainty surrounding seasonal occurrence, overlapping O&G platforms, and primary foraging technique were all lower, at 0.2. Flocking behavior had the lowest uncertainty (0.10) ( Table 2 ).

The number of detections for each species ranged from 27 - 1,104, and the AUC for Maxent models ranged from 0.872 – 0.996 (training dataset) and 0.815 – 0.975 (testing dataset), suggesting very good to excellent model performance (Swets, 1988; Duan et al., 2014 ; Table 4 ). The spatial extent of highly suitable habitat ranged from ~9,100 to ~281,700 km² of the study area among the focal species ( Table 4 ; Figures S3A–X ).

The percent of species-specific habitat within the O&G platform footprint varied among species (range = 0% - 64.96%, median = 14.00%, mean = 21.67% ± 13.97% SE; Table 2 ). The overlap of highly suitable habitat with the O&G platform footprint was greatest for Brown pelican (Pelecanus occidentalis; 64.96%) and least for Brown noddy (0%) ( Table 2 ; Figures S3 I, H ). Species with mostly nearshore habitat, such as Brown pelican and Black tern, had the highest percent of habitat within the O&G platform footprint ( Table 2 ; Figures S3I, C ). Species with highly pelagic distributions, such as Brown booby and Cory’s shearwater (Calonectris borealis), had a lower percent of habitat O&G platform footprint ( Table 2 ; Figures S3G, l ). The lowest percent of habitat within the O&G platform footprint occurred in species with most of their habitat in the eastern portion of the study area, such as Sooty tern and Brown noddy ( Table 2 , Figures S3W, H ).

The seabird habitat footprint covered ~521,600 km², encompassing portions of all bathymetric domains and covering much of the study area ( Figure 4 ). The total area of the O&G platform footprint (platform locations surrounded by a 10 km buffer) was ~88,400 km². Much of the O&G platform footprint occurred in waters < 500 m depth. Approximately 15% of the seabird habitat footprint (~78,900 km²) overlapped the O&G platform footprint, mainly < 500 m depth in the central and, to a lesser extent, western portions of the study area. A relatively large portion of the O&G platform footprint, ~89%, overlapped the seabird habitat footprint. This suggests that most O&G platforms co-occur with the habitats of one or more seabird species. Due to the absence of O&G platforms, no overlap occurred in the eastern portion of the study area ( Figure 4 ).

Cumulative seabird vulnerability was high along the continental shelf-slope in waters between ~200 m and 2,000 m ( Figure 5 ). Moderate to low cumulative species vulnerability occasionally extended into pelagic areas. The cumulative seabird vulnerability/km² was 3.21 within the O&G footprint compared to 2.33 outside of the O&G footprint. Overlap with high and moderate cumulative seabird vulnerability occurred from nearshore areas to the mid-shelf where the depth was ~500 m, and into more pelagic waters south of Louisiana.

The two most vulnerable species, Northern gannet and Audubon’s shearwater, had similar life-history and behavioral variable scores but very different levels of spatial overlap with O&G platforms (i.e., Northern gannet ~58%, Audubon’s shearwater ~5%). Breeding in Labrador and Newfoundland, ~25% of Northern gannets migrate to the nGoM during the nonbreeding/wintering period (November) and generally remain in the Gulf until they depart in the spring (April) (Montevecchi et al., 2012; Fifield et al., 2014). Blood, feather, and isotopic analysis of Northern gannet wintering in the nGoM indicate that they experience more stress than birds originating from the same colony, Bonaventure Island (Quebec, Canada), but wintering on the U.S. Atlantic coast (Champoux et al., 2020). The drivers of these differences are uncertain, but they may relate to greater exposure to environmental stressors in the nGoM than in the U.S. South Atlantic coast. Although Northern gannet are considered a species of least concern globally, our data and others suggest a relatively high level of vulnerability in the northern Gulf of Mexico (Jodice et al., 2019).

In contrast to Northern gannet using the nGoM in the nonbreeding season, the breeding origin of Audubon’s shearwater using the nGoM is not well understood. Audubon’s shearwater colonies exist throughout the Caribbean, but birds tagged in colonies in the northern Bahamas and Martinique did not enter the nGoM (Precheur, 2015; Ramos et al., 2021). Cay Sal Bank, which supports ~5,000 breeding pairs, is the closest known breeding site to the Gulf of Mexico (Mackin, 2016). With the most detections occurring in spring, summer, and fall, it is likely that many of the Audubon’s shearwater observed in the nGoM were post-breeding individuals. In addition to the variables included in the index applied here, the use of Sargassum reefs by Audubon’s shearwater for foraging (Haney, 1986; Moser and Lee, 2012) could lead to an additional oil exposure pathway, as algal mats and oil can aggregate due to processes at the water’s surface (Carmichael et al., 2012; Powers et al., 2013). While globally considered a species of least concern, Audubon’s shearwater has a Partners in Flight (United States of America and Canada) watch list classification of Yellow-D, indicating steep declines and the existence of major threats, and the Caribbean subspecies likely using the nGoM is a Caribbean at-risk species (Bradley and Norton, 2009; BirdLife International, 2021; Partners in Flight, 2021). A better understanding of the origins and population trends for Audubon’s shearwater using the nGoM could inform potential research or restoration actions related to this species.

The high uncertainty of sensitivity variables is not surprising as the data needed to inform these variables can require intensive, multi-annual observation and monitoring efforts, and this is often lacking for seabirds that nest in remote locations. For example, if the determination of age at first breeding was deemed a variable for which additional data were required to fill a data gap, then a sustained population monitoring effort would likely need to be employed across a wide geographic range. Likewise, nest-based monitoring could provide additional data on the number of days from incubation through fledge. Beyond reducing the uncertainty in such variables, improved estimates of demographic data would allow other metrics synthesizing multiple life-history characteristics to be estimated. Such data could estimate the number of years lost through adult mortalities and the subsequent loss of potential offspring over their expected lifetimes, i.e., ‘bird years’ (Cole and Dahl, 2013; Evers et al., 2019). Within a Resource Equivalency Analysis framework, lost bird years have been used as a currency to estimate the degree of compensation needed to offset bird mortality associated with oil spills (Evers et al., 2019). Therefore, collecting robust demographic information could reduce uncertainty in sensitivity variables and augment the available tools for estimating interactions’ impacts and defining compensation targets.

Uncertainty in residency could be reduced by increasing the understanding of seabird movements across their annual cycle, including the degree of individual and interannual variation in movement patterns. Lagrangian methods, such as individual-based tracking, can be used to better understand variation in migration patterns (Baak et al., 2021), habitat use (Jodice et al., 2015; McCloskey et al., 2018; Phillips et al., 2018), and overlap with threats (Lamb et al., 2018; Isaksson et al., 2021). Tracking has identified previously unknown breeding areas (Kanai et al., 2002) and demonstrated links with distant colonies with seabirds using the nGoM (Montevecchi et al., 2012). Such information could inform population monitoring or recovery efforts and direct them towards the colonies and populations known to use the nGoM. More information on the magnitude, frequency, and degree of overlap with threats outside of the nGoM (e.g., bycatch, threats at breeding colonies), would enhance the inference gained from this variable, potentially revealing unexplored options to address seabird interactions in the nGoM. As data from the GoMMAPPS program are assessed, details on the marine range of species occurring in the Gulf are emerging (e.g., Jodice et al., 2021). Given that relatively few vessel surveys and satellite tracking efforts exist in the nGoM compared to other U.S. oil and gas production regions (e.g., Atlantic coast, California, Alaska), sustained, long-term regional surveys and satellite tracking would reduce uncertainty in estimates of seabird vulnerability to oil.

The extensive seabird habitat footprint demonstrates broad-scale use of the nGoM by seabirds, and the disproportionate percent of cumulative seabird vulnerability overlapping O&G platforms, notably on the shelf-slope, illustrates the non-exclusive use of Gulf habitats by humans and wildlife. Shelf-slope areas may offer seabirds enhanced foraging opportunities and prey aggregation associated with frontal zones such as the Louisiana Texas Shelf Front or other mesoscale features. As convergence zones, ocean fronts, and other dynamic features have been associated with seabird foraging behavior (Haney and McGillivary, 1985; Weimerskirch, 2007; Scales et al., 2014; Poli et al., 2017), these features could have a substantial impact on seabird distribution. Further study of seabird associations with oceanographic features in the nGoM could investigate these potential relationships and refine the characterization of seabird habitat, which could be used to reduce, manage, or mitigate interactions with O&G activities.

The near-complete (~89%) overlap of O&G platforms within the seabird habitat footprint suggests a high potential for seabirds to interact with any given platform. Much of this overlap of seabird habitat with O&G platforms in waters< 500 m depth, particularly near the 200 m isobath, an area of high cumulative seabird vulnerability. Seabird surveys in the Gulf of Mexico have noted relatively low densities of seabirds over the middle and outer continental shelf (Haney et al., 2019). However, without detailed information on species-specific densities, it is unclear how the number of individual seabirds could impact the characteristics of the vulnerability of the entire assemblage of seabirds in the nGoM.

In addition to small spills and produced waters, seabirds and other avifauna can be impacted by platforms themselves. Platforms and platform activities, particularly lights and flares, can distract and disorient avifauna, potentially resulting in collisions (Russell, 2005; Montevecchi, 2006; Ronconi et al., 2015). A better understanding of the spatial extent of seabird attraction to flares, which can occur at night and during the day, could be used to refine how the spatial footprint of O&G is defined and improve our understanding of the impacts of O&G activities. Seabirds in the nGoM have been observed roosting and foraging from platforms (Ortego, 1978; Russell, 2005), potentially resulting in prolonged/increased exposure to oiling and other anthropogenic interactions. Direct and indirect interactions with O&G platforms could be better understood through standardized monitoring programs with trained observers, systematically covering different environmental and lighting conditions (Wiese et al., 2001; Burke et al., 2012; Ronconi et al., 2015). Technologies including telemetry, cameras, and radar could also be used to augment monitoring efforts and coverage (Ronconi et al., 2015). While overlap at a macro-scale does not equate to an interaction, the degree of O&G platforms overlapping seabird habitat in our study area merits consideration as a coarse indicator of potential interaction.

Ongoing development of O&G activity in the nGoM suggests that seabird overlap and potential interactions with offshore energy activities are not likely to decline over time. Continued installation of ultra-deep wells (> 1,000 meters depth) will expand the O&G footprint along the continental shelf and into more pelagic waters (Murawski et al., 2020). This spatial expansion would increase the spatial extent of overlap with seabird habitats. Increased overlap would increase seabird vulnerability scores of many seabird species, particularly those that occupy pelagic habitats. An expanded footprint of O&G platforms would also augment acute oil exposure via spills and chronic exposure through the continued release of produced waters. Thus, updating the scores of variables in the index to reflect current levels of overlap could help maintain the relevancy of index-based relative vulnerability ranks. Some have suggested that climate change could amplify these risks, with increased storm frequency heightening the potential for damage to oil and gas infrastructure and spills (Cruz and Krausmann, 2013). However, the relative importance of mechanisms controlling tropical cyclone activity in the Gulf remains unclear; the direction of change (increase or decrease) of future tropical cyclone activity in the nGoM remains an active field of research (Colbert et al., 2013; Knutson et al., 2015; Rodysill et al., 2020).

Installation of wind turbines on the Gulf of Mexico Outer Continental Shelf (DOI, 2021) could increase the footprint and density of marine infrastructure in areas ≤ 500 m depth, the current maximum depth for wind turbines. This depth range coincides with areas of high cumulative seabird vulnerability. The behavioral response and impacts of offshore wind energy development will be species-specific and may involve avoidance or attraction (Furness et al., 2013; Dierschke et al., 2016). Avoidance behavior can increase energy expenditure, while attraction can lead to collision mortalities or other lethal or sub-lethal interactions (Drewitt and Langston, 2006; Masden et al., 2010). These behavioral modifications and impacts could be compounded with those related to O&G activities. Understanding the cumulative effects of offshore energy, characterizing how and when seabirds interact with offshore energy acquisition, and considering the impacts of climate change appears warranted.