An ecological network approach to predict ecosystem service vulnerability to species losses

Human-driven threats are changing biodiversity, impacting ecosystem services. The loss of one species can trigger secondary extinctions of additional species, because species interact–yet the consequences of these secondary extinctions for services remain underexplored. Herein, we compare robustness of food webs and the ecosystem services (hereafter ‘services’) they provide; and investigate factors determining service responses to secondary extinctions. Simulating twelve extinction scenarios for estuarine food webs with seven services, we find that food web and service robustness are highly correlated, but that robustness varies across services depending on their trophic level and redundancy. Further, we find that species providing services do not play a critical role in stabilizing food webs – whereas species playing supporting roles in services through interactions are critical to the robustness of both food webs and services. Together, our results reveal indirect risks to services through secondary species losses and predictable differences in vulnerability across services.


I.1 Ecosystem Service Assignment to Species
We assigned ecosystem services to species after filtering the original salt marsh food webs 1 to exclude parasites and detritus, and to only include adult life stages. Then, we conducted a literature search in Google Scholar to select the main ecosystem services provided by salt marsh ecosystems: shoreline protection, wave attenuation, water filtration, carbon sequestration, fisheries, birdwatching and waterfowl hunting. To determine ecosystem service providers for shoreline protection, wave attenuation, water filtration, carbon sequestration, we conducted a literature search in Google Scholar, following previous methods 2 . Specifically, we searched for each unique servicespecies combination, using both species common and scientific names, as well as synonyms for each ecosystem service. For example, carbon sequestration could also be published as "carbon storage." If a reference existed that indicated a species provided a service, we augmented the interaction data (link data) to include a link from the species to the service. For the fishery and waterfowl ecosystem services we referenced government documents with seasonal regulations and/or catch limits and guidelines (e.g. 3,4 ). Data will be published for the ecosystem service network layer (Dee and Keyes, in review).
Birdwatching. We used the R package auk v 0.4.1 to process and analyze the data 5 . We filtered the eBird Basic Dataset by Complete Checklists and by spatial extent matching the eBird registered hotspots with the latitude and longitude from the original food web data 1 . We did not filter by month or year, instead using all records from 1900-2019. For each estuary, we aggregated count data by species to get a total count (i.e., total number of birds seen per species across all complete checklists) which we used to calculate the relative frequency at each estuary. Species with lower relative frequencies were referred to as rarer than those with higher relative frequencies.
Waterfowl Hunting. For Carpinteria Salt Marsh, we assigned links for bird species that are hunted and regulated in the Southern California Zone, which includes Santa Barbara, USA. For Bahia San Quintin and Estero de Punta Banda, MX, we assigned links for birds listed with hunting seasons and limits in the SEMARNAT working season document 6 .
Fisheries. For the Carpinteria Salt Marsh, we referenced California Saltwater Sport and Angling fishing regulatory documents 3,7 . For the two salt marshes in Baja, Mexico, we referenced the National Fishing Charter provided by the Instituto Nacional de Pesca 4 . We assigned fishery links to any species that was reported as fished in the northern region of Baja, Mexico.

I.2 Species Loss Sequences
Ecosystem Service Providers: The species removed in this extinction sequence are the species directly providing ecosystem services. These include species that were identified in section 2.2 that were linked to services directly in the link data. For example, birds that were linked to waterfowl hunting and fish that were linked to fishing because they were listed in a regulatory document were removed in this extinction sequence. We first ran this sequence three times using different orders of species loss before deciding to use the biomass-ordered sequence discussed in the main text. We first filtered the interaction data to only include the species that provide services directly. We further simplified the data so that species that provide multiple services were only listed once. We simulated the loss of direct ecosystem service providers in three ways: 1. By biomass, in ascending order, 2. By abundance, in ascending order, and 3. Randomly, simulated 100 times. We estimated missing abundance and biomass data using previously published formulas 8 . Following this, we had biomass data for all species and abundance data for most species. Both food web and ecosystem service robustness were higher under the abundance-ordered sequence (see Figure S4). We chose to include the biomass-and random-ordered sequence in our final results because we had more data for these sequences. This extinction sequence did not include all the species, so the simulations ended when the last service provider was removed.

Figure 2.
Robustness values vary for the removal of ecosystem service providers by biomass (low to high biomass) and abundance (least to most abundant) for A) the aggregate ES robustness (RES) and B) food web robustness (RF). Results are shown for each system: BSQ, CSM, and EPM.

Supporting species:
We calculated the personalized PageRank 9,10 score for each species in the food web, once for each ecosystem service ( Figure S3 and S4). An important specification in a PageRank algorithm is αthe probability that the random walker terminates its walk at any given node ( Figure S4). The magnitude of α determines how far the random walker gets into the network. The higher the value of α, the more likely the random walker terminates its walk closer to the starting node. Commonly, α is set to 0.15 9 . We tested various values of α to assess potential variation in the "importance" of supporting species, i.e. species with a high probability of the walker visiting that node. We calculated the personalized PageRank of all species for each ecosystem service using the following α values: 0.15, 0.35, and 0.45. We then calculated each species mean contribution to ecosystem services in aggregate (see main Methods for details) and calculated the robustness of the food webs and ecosystem services with the following sequences of species losses: high to low importance and low to high importance (for all three α values). We did not observe substantial variation in food web or ecosystem service robustness values (except see BSQ Food Web, α = 0.5; Figure S3), thus we followed commonly used α values and included the results from α = 0.15. Personalized PageRank to identify critical supporting species. Arrows point from ecosystem services to ecosystem service providers, and from consumers to resources. (B) Personalized PageRank allows us to specify the starting node as ecosystem service nodes. The walker starts every walk at the ecosystem service node, pursuing paths through the network at random. At each new node, the probability that the random walk terminates at a given node is α, and the probability that the walk continues is 1-α. Here, in this simplified diagram, 2/3 of the random walks end at node b, and 1/3 of the walks end at node a, indicating that node b would have a higher probability than a, and both have a higher probability than the rest of the nodes. The higher probability can be interpreted as "more important" to the ecosystem service node.  . Alpha is the probability that the random walk terminates at any given node (α in Figure S3). Alpha (α) values ranged from 0.15-0.5, where values closer to 1 are high and values closer to 0 are low. The higher the value of alpha, the shorter the random walks (i.e., the random walker may not reach the basal trophic levels as frequently).

Rarity and body size:
Previous food web robustness studies have simulated species loss ordered based on body size [11][12][13] , which can influence extinction risk 14 . However, both large-bodied species and smallbodied species can be vulnerable to disturbances. We ran two body size sequences following the literature (removing from large to small and small to large bodied species 14,15 ). We used body size information from the original food web data 1 . However, this sequence did not include plants because there is not plant body size information. To account for this, we instead included the Rarity sequence, where species are removed from low to high relative abundance (estimates based on Hechinger et al 1, 8 ). The order of the sequence was identical except for the difference in inclusion of plant species, so we decided to only include the rarity sequence in the main text.

Species impacted by system-specific threats, 'vulnerable species':
Prior to conducting the literature synthesis to determine which species are threatened, we performed a literature search to identify salt marsh-specific threats. We searched Google Scholar using the following terms: salt marsh, estuary, disturbance(s), threat(s), risk(s). Once we established estuarine acidification, eutrophication and hypoxia, and pollution as our disturbance we searched for species in our food webs that are known to be impacted by them. We looked at each species individually using a combination of its scientific or common name and each of the disturbances. If a reference existed to support that species s displayed known mortality to the threat i, "1" was recorded to indicate the species was threatened, and the reference was recorded. We recorded "0" for all species that did not have a record of known mortality. While this method allowed us to systematically identify threatened species, we likely underestimated the number of threatened species. Some species are part of taxonomic groups that have varying impacts to these disturbances, such as fish.
For these species, we only recorded "1" if the specific species name displayed known mortality. The species threatened by these stressors was a small subset of the full species list. Supplementary Data 1 includes the threatened species across the three estuaries and the reference supporting its categorization.

I.3 Robustness methods and sensitivity analyses
We simulated extinctions in R, relying on the package purr in tidyverse v 0.3.4 16 . The procedure is as follows: a target species is removed from the network, then any species who no longer has any resources in the network is considered to be secondarily extinct and are tallied. The process of target species removal continues until there are no remaining target species in the sequence.
In the main text, we calculate robustness as the area under the curve where the x-axis is the proportion of target species removed and the y-axis is proportion of species/ecosystem services that did not go secondarily extinct. We calculated the y-axis as: − where S is the number of nodes susceptible to secondary extinction (i.e., non-basal) and C is the cumulative number of secondary extinctions at each time step, t. Species nodes are susceptible if they have at least one resource (in-degree > 0, i.e., non-basal species), and all ecosystem services are susceptible because they inherently rely on at least one species (all in-degree > 0). This calculation allows us to only track secondary losses, regardless of whether the node represents a species or an ecosystem service. If the y-axis does not reach 0 in the food web, that indicates that some susceptible species were directly removed (i.e., primary extinctions on the x-axis), and/or that there are susceptible species remaining in the food web following the completion of a sequence.
We tested the sensitivity of our results to the choice of x-axis for calculating the Area Under the Curve (following 17 ). Specifically, we examined multiple variations for the x-axis to make sure that the correlation between food web and ecosystem service robustness was consistent. The robustness values (AUC) varied based on the x-axis calculation, so it was important to test that the relationship between the two values was not sensitive to this calculation.
We calculated the robustness of the Carpinteria Salt Marsh (CSM) food web and ecosystem services for all 12 sequences using the following two equations to calculate the x-axis: (1) where D is the number of species directly removed (i.e., number of realized primary extinctions), and LSeq is the length of the intended sequence; and (2) where D is the number of species directly removed (i.e., the number of realized primary extinctions), and N is the total number of species in the food web.
We then ran two Spearman's correlation tests on the food web and ecosystem service robustness for CSM, once using robustness values calculated with x-axis equation (1) above, and once using robustness values calculated from x-axis equation (2).
The main results remained unchanged: Food web and ecosystem service robustness are strongly and positively correlated for both equation (1) (ρ [12]=0.93, p < 0.001) and equation (2) (ρ [12]=0.881, p < 0.001). We used robustness values that used equation (1) for the x-axis because it better reflects the response of the food webs and ecosystem services to the specific order of species removals in each of the sequences.  and (c), the x-axis is proportion of target species removed, which increases towards 1 every time a species is removed in a species loss sequence. In (a), the y-axis is proportion of species remaining following the removal of a species. In (b), the y-axis is the proportion of services (in aggregate, including all 7 services) remaining following the removal of a species. Each decrease in the y-axis shows when an ecosystem service was lost. In (c), the y-axis indicates whether a single ecosystem service remains (y=1) or was secondarily lost (y=0). The key difference is that in food web robustness, we are tracking secondary species losses that result from primary species removal, and in ecosystem service robustness, we are tracking ecosystem service losses (either in aggregate (b) or singularly (c)) that result from primary species removals.

Robustness of ecosystem services with weighted species contributions
Species contributions to many ecosystem services scale proportionally to their biomass (e.g., regulating and supporting services [18][19][20] , but see 21 ). For instance, species with higher biomass contribute more to carbon fixation and nutrient cycling than species with low biomass 22 . We extended the ecosystem service robustness calculations in the Robustness Analysis section (main text) to account for differences in species' contributions to each service, proportionate to an ecosystem service provider's biomass. First, we estimated biomass for species without data because the original data 1 do not include biomass estimates for all species. To do this, we used methods explained in 8 where biomass is a function of species' abundances and body sizes: Then, we simulated the same 12 extinction scenarios but tracked the proportional change in ecosystem services ( Figure S6(a)) as opposed to the binary outcomes of retention or loss (i.e., complete service loss when all ecosystem service providers are lost). We calculated proportional changes in ecosystem service provision as an aggregate ( Figure S6(a)) and as individual services ( Figure S6(b)). The size of all ecosystem services was considered to be equivalent to the summed biomasses of all the species directly providing them. The key difference between the aggregate and individual weighted robustness is that the aggregate pools all ecosystem service providers' biomasses together, regardless of the service; and the individual calculation only considers species that are directly linked to a single service. For the aggregate, the proportion of ecosystem services remaining at each time step was calculated using: With this calculation, each time an ecosystem service provider is lost, the y-axis decreases proportionally to its biomass. We calculated the weighted, individual ecosystem service robustness in the same way, except the y-axis only considers the ecosystem service providers directly linked to that service ( Figure  S6(b)). So, the denominator is the total ecosystem service provider biomass for the single ecosystem service, and the numerator is the remaining ecosystem service biomass for that single service (i.e., 1cumulative sum of ecosystem service provider biomass removed or secondarily lost). Each time an ecosystem service provider is lost, the y-axis decreases proportionally to its biomass.
We did not do weighted contributions of species for birdwatching because the value of this service does not scale proportionally with a species' biomass 19,21,23 .

Figure 6 Illustration of weighted robustness calculations for (a) ecosystem services in aggregate and (b) individual ecosystem services. In (a) and (b)
, the x-axis is the proportion of target species removed, which increases towards 1 every time a species is removed in a species loss sequence (as in the unweighted robustness calculation in X). In (a), the y-axis is the proportion of all services in aggregate (including all services except birdwatching) remaining following the removal and/or the secondary loss of an ecosystem service provider. Each decrease in the y-axis shows the percent of all ecosystem services that were lost, which is based on the biomass of the ecosystem service providers that were removed or secondarily lost. In (b), the y-axis indicates decreases in a single ecosystem service, which is based on the proportion of all ecosystem service providers biomass for that single ecosystem service. Each time an ecosystem service provider is lost for that ecosystem service, the y-axis decreases proportionally with that species' biomass.   Figure 1A), and Estero de Punta Banda ( Figure 1B). The color denotes the sequence of species removals: ecosystem service providers (removed by biomass and randomly), supporting species (removed by high-to-low importance, low-to-high importance, and randomly), most-to-least connected, randomly (reporting: highest, lowest and mean), rarity (removed by relative abundance -from least to most abundant), and vulnerable species (removed from least to most abundant).

II.1 Food web and ecosystem service robustness strong, positive correlation not skewed by the length of the vulnerable species sequence.
We ran sequences (i.e., species loss scenarios) that varied in length, where some ended after very few species were removed and others ended when all species were removed. The vulnerable species sequence was short (Lseq ranged from 5-20 across the 3 salt marsh systems), due to our conservative approach to include species in the sequence (see main text). This short sequence length could mirror reality because it is not realistic that all species are always lost. The removal of vulnerable species (based on our methods for assigning threat status) did not cause many (if any) secondary extinctions, resulting in high food web and ecosystem service robustness. However, we wanted to be sure that these high (and very similar) robustness values did not skew our correlation tests (e.g., artificial strong correlation). To test this, we ran two correlation tests on the food web and ecosystem service robustness values from Carpinteria Salt Marsh (CSM): (1) data from all sequences, including the vulnerable species data, and (2) data from all sequences, excluding the vulnerable species data. The main results remain unchanged. The correlation between food web and ecosystem service robustness was strong, positive and significant (with p < 0.001) when the vulnerable species sequence was included (ρ[36]=0.884) and excluded (ρ[33]=0.849), so we included it in our main results.

II.2 Ecosystem service robustness across sequences and systems
Ecosystem service robustness was varied across all sequences (min. RES = 0.256, max. RES = 1; main Figure 4(a)). Ecosystem service robustness was highest (RES = 1) for all three systems under the systemspecific threats sequence because no services were lost after all species were removed. The systemspecific threats sequence was the shortest sequence (with 5-20 species directly removed). Ecosystem service robustness was lowest for all three systems under the most-to-least important supporting species sequence (min. RES=0.256, max. RES=0.320). Ecosystem service robustness varied across systems (difference in min. and max. RES ranged from 0-0.261).

II.3 Food web robustness across sequences
Food web robustness varied across all 12 sequences of species loss (min. RF = 0.145, max. RF = 1; main Figure 4(b)). Under the most-to-least connected and the most-to-least important supporting species sequences, the food web collapsed before more than 50% of the species were removed directly (main Figure 5C) Similarly, the random removal of species resulted in secondary species losses sooner than the threat-based sequences and the important ecosystem service provider sequence.

II.4 Sensitivity of results to sequences that removed different numbers of species
We simulated sequences of varying lengths, ranging from 5 species to all species. Table S2 shows whether or not a sequence included the full length of the species list. We only compared the raw robustness value for sequences that included the full species list. To compare the short sequence, 'Vulnerable species', to the other sequences of species lost, we compared the raw number of species that were secondarily lost in the removal of the same number of target species for the other sequences. We found that there were varying levels of secondary species losses across the sequences in the same number of target removals (Table S3).

II.5 Is ecosystem service robustness a function of sample size?
Given that ecosystem service robustness differed from food web robustness, and that this approach has not yet been applied, we explored the extent to which our ecosystem service robustness results was a function of sample size (three salt marsh systems, 7 ecosystem services each). By matching each ecosystem service to a similar species (in-degree centrality and trophic level), we are able to explore how are results compare to a common unit in robustness studiesspecies.
We calculated the in-degree centrality and trophic level of all nodes in the CSM food web-ES network using iGraph v 1.2.5 10 and NetIndicies v 1.4.4 24 R packages, respectively. Once we had this information, we identified 2 species nodes that closely matched the in-degree centrality and/or the trophic level of each of the 7 ecosystem services (Table S4). We then ran the same 14 robustness sequences (but with these 7 matching species removed from the sequences) and calculated the overall robustness of the 7 identified species (in aggregate, Table S5). Finally, we calculated the mean number of species remaining when each ecosystem service was lost, and the mean number of species that were lost prior to losing each ecosystem service. We did this for the 7 matched species as well (Table S6).
From this sensitivity analysis, we see that species similar to ecosystem services (in terms of in-degree centrality and trophic level) behave similarly to ecosystem services themselves when species are removed in different orders. However, the robustness of species similar to ecosystem services is higher or lower than ecosystem service robustness depending on the sequence of species removals. Further, ecosystem service nodes and the species similar to ecosystem services were lost at different times in each sequence (Table S6, Figure S8). They similarly display threshold behavior ( Figure S8).   Table 6. Tracking when ecosystem service nodes and matching species nodes are lost. For each ecosystem service node and matching species nodes (target nodes), we tracked how many species were lost and remaining when each target node was secondarily lost. Using these raw numbers for each of the 12 sequences, we calculated the mean proportion of species lost and remaining, where the denominator is the total number of species in the system. This allowed us to compare the sensitivity of the ecosystem service nodes and the matching species to each sequence of species loss.  Figure 10. Species losses causing the secondary loss of ecosystem services and species similar to ecosystem services. We calculated the robustness of ecosystem service nodes and species nodes with similar trophic level and in-degree centrality using the data from these plots. Each box shows the sequential loss of species (x-axis) and the secondary losses of ecosystem services (y-axis, blue) or species similar to ecosystem services (y-axis, yellow).

II.6 Sensitivity of results to the exclusion of parasites
Parasites play an important and unique role in estuarine foods webs 25,26 , altering their structure 27 and stability 28 . Previous studies in the three estuaries under consideration here 1 have demonstrated that the inclusion of parasites can decrease network robustness 29 . In these estuaries, the susceptibility of parasites to secondary extinctions is a result of their complex life histories which require them to pass through a series of hosts in specific sequence in order to complete their life cycle. Parasites in these estuaries are dominated by trematodes, many of which must pass through a single species of snail as larva, but which can occur in dozens of species of fish-eating birds as adults. Thus, the parasite influences on robustness are a result of ontogenetic diet shifts (often serial specialization) which can be captured by stagestructured analysis. While parasites have stage-specific diet information in our networks, their free-living counterparts (which comprise the majority of species) do not. This precluded stage-structured analyses. Given their potential robustness impacts, we incorporated adult parasites and re-ran our robustness analyses for the food webs and ecosystem services for all sequences on networks that included parasites and parasitic interactions (see Table S8). To examine the sensitivity our results to including or excluding parasitic interactions, we ran a Spearman's correlation test on the food web and ecosystem service for all 3 systems and all 12 sequences.  Figure S11). Generally (for more than 6 of the 12 sequences), robustness values were higher when parasites were included for both the food webs and the ecosystem services (across all three systems; Figure S12). This is due to the fact that we only included adult life stages, which for these parasites tend towards diet generality inflating robustness 25,29 . Interestingly, food web robustness was substantially higher when parasites were included in all three salt marshes when supporting species were removed from most to least important ( Figure S12: 'I' Food Web). However, similar to previous studies that have compared food web robustness with and without parasites on these food webs, robustness decreased for both food webs and ecosystem services when parasites were included for two sequences: most to least connected 26 and random 28,29 . Figure 11. Food web and ecosystem service robustness remain highly and positively correlated when parasites are included (rs[36] = 0.738, p < 0.001). A comparison of food web vs. ecosystem service robustness across each species loss scenario and salt marsh system. The black line shows the linear regression line. The shape of data points represents the salt marsh system: Bahia Falsa de San Quintin, Carpinteria Salt Marsh, and Estero de Punta Banda. The color denotes the sequence of species removals: ecosystem service providers (removed by biomass and randomly), supporting species (removed by highto-low importance, low-to-high importance, and randomly), most-to-least connected, randomly (reporting: highest, lowest and mean), rarity (removed by relative abundance -from least to most abundant), and vulnerable species (removed from least to most abundant).  Table 8. Coefficient estimates and model summary for each of the 7 linear models. We fit 1 model of each sequence B-D and 1 model that included the data for all sequences. We did not make adjustments for multiple comparisons. The dependent variable in each model was individual ecosystem service robustness (Rindiv) and the independent variables were mean trophic level and link redundancy. Model: A) All sequences, B) Most-to-least connected, C) Ecosystem service providers (high-to-low biomass), D) Ecosystem service providers (low-to-high biomass), E) Supporting species (most-to-least important), F) Supporting species (least-to-most important), G) Rarity.