Co-extinction in a host-parasite network: identifying key hosts for network stability

Parasites comprise a substantial portion of total biodiversity. Ultimately, this means that host extinction could result in many secondary extinctions of obligate parasites and potentially alter host-parasite network structure. Here, we examined a highly resolved fish-parasite network to determine key hosts responsible for maintaining parasite diversity and network structure (quantified here as nestedness and modularity). We evaluated four possible host extinction orders and compared the resulting co-extinction dynamics to random extinction simulations; including host removal based on estimated extinction risk, parasite species richness and host level contributions to nestedness and modularity. We found that all extinction orders, except the one based on realistic extinction risk, resulted in faster declines in parasite diversity and network structure relative to random biodiversity loss. Further, we determined species-level contributions to network structure were best predicted by parasite species richness and host family. Taken together, we demonstrate that a small proportion of hosts contribute substantially to network structure and that removal of these hosts results in rapid declines in parasite diversity and network structure. As network stability can potentially be inferred through measures of network structure, our findings may provide insight into species traits that confer stability.


Materials and Methods
Study system. The current study utilizes data from a Long Term Ecological Research study on the fish and parasite communities of the Upper Paraná River floodplain, a 230 km long dam-free portion of the Paraná River in southern Brazil. Fish were sampled quarterly for eight years (2000-2008) by using both active and passive sampling techniques, resulting in a total of 4,875 fish captures of 72 species selected for parasitological analysis, Native species (n = 49) outnumbered non-native species (n = 23) of those fish species chosen for parasitological analyses. Parasites (n = 324) were identified to the best possible resolution, which resulted in 140 parasites identified to the species level, 120 to genus, and 64 to family 24 . The network is significantly nested and modular 24 , and has low connectance (C = 0.022), similar to other fish-parasite networks 25 . Details for the sampling protocols 26 and parasite examination 27,28 are provided elsewhere. Calculations of network structure metrics were based on the host-parasite interaction matrix provided by Lima et al. (2012) 24 .
Fish host trait variables selected. Data on fish traits (Table 1) were obtained from FishBase [Froese, R. and D. Pauly, www.FishBase.org, accessed 08/2014] a comprehensive fish life history trait database accessible in R using rfishbase 29 , and based on information provided in Lima et al. (2012) 24 . The importance of certain fish species in the network may relate to specific host traits. For instance, body size and range size are often correlated, and larger ranged hosts typically harbor higher parasite diversity. Loss of these hosts could have a larger negative impact on network structure. However, we did not exclude variables for which we had data on, as there was no a priori reason to exclude data, and our analysis method was insensitive to issues of collinearity and overparameterization. Parasite species richness was quantified as the number of parasite species found on a given host species. Measures obtained from FishBase included host extinction risk ("vulnerability") 30 , resilience from a perturbation 31 , phylogenetic diversity 32 , maximum weight, body size, length at maturity, and host family. Vulnerability is a predicted measure of extinction risk based on host life history information. Resilience is a measure that incorporates host demographic information (carrying capacity, fecundity, etc.) to classify the likelihood that a fish population will recover quickly from a perturbation. Phylogenetic diversity is reported in FishBase and based on Faith et al. (2004) 32 . The specific measure used, PD 50 , represents the expected loss in phylogenetic diversity under species extinction, assuming that all other species have equiprobable extinction risk. Measure of sampling effort, fish abundance, and biomass were obtained from Lima et al. 2012 24 . Host species "status" (i.e., native or non-native) was compiled by Lima et al. (2012) 24 , though the original evaluation was performed by Julio et al. (2009) 33 .
Host species-level contribution to nestedness and modularity. We quantified nestedness using a metric based on overlap and decreasing matrix fill (NODF) 34 . Nestedness contribution (C i ), a measure of the importance of a single node to overall network structure, was calculated following Saavedra Scientific RepoRts | 5:13185 | DOi: 10.1038/srep13185 (2011) 18 , in which the interactions of a single host species were iteratively randomized, with the resulting simulated NODF values compared to the observed nestedness of the empirical network. Modularity contribution (Q i ) was calculated in the same way, but using Barber's Q statistic 20 . Both measures were calculated only for the complete network. It is likely that species contribution to nestedness and modularity changes as the network disassembles. However, identifying nestedness and modularity contributions via dynamic updating makes the measures state-dependent, in such that the relative impact of a species removal in the full network may be greatly altered after a certain number of nodes have been removed. As a result, our modularity contribution values do not contain information on how the network will respond when host species are lost from the network. Values of C i or Q i close to zero indicate little or null impact of that host on network structure, while positive (negative) values are indicative of hosts that enhance (reduce) nestedness or modularity.
We used three randomization algorithms ranging from a) each interaction equiprobable b) probability of interaction determined by parasite host breadth and c) probability of interaction determined by average of host and parasite occurrence probabilities (supplementary equation S1; Saavedra 2011 18 ). All null models gave similar results, so we report parasite host breadth to determine occurrence probabilities herein (b from above), and relationships between measures obtained from the different null models are presented as supplementary material (Figs S1 & S2).
Traits associated with host nestedness contribution. In order to assess the potential relationship between host traits and their corresponding contributions to network structure (i.e. C i and Q i ), we used regression tree analysis. Regression trees are a useful tool for the identification of the most important predictor variables, and have been used extensively in ecological applications [35][36][37][38] . Regression tree analysis aims to predict a response variable based on several input variables, in which input variables represent splits in the tree. We applied an iterative process called "boosting", where we formed many regression trees, and then combined to improve predictive power. The final boosted regression tree (BRT) model was cross-validated using ten-fold cross validation and permutation tests to determine the optimal number of trees. Variables included in the BRT analysis are provided in Table 1. The relative contribution (RC) of predictor variables was determined through permutation tests, in which each predictor variable was permuted randomly, and the importance of the predictor variable was quantified by the subsequent reduction in model performance. Relative contribution estimates are based on the number of times a given predictor variable is selected for splitting, weighted by the degree the split improves model performance. This metric is then averaged across all trees built in the model, and scaled between 0 (no contribution) to 100 (high contribution).
Host removal simulations. In order to assess the impact of host extinction on network structure and parasite diversity, we sequentially removed hosts according to four potential extinction orders. We removed hosts from the network in order from highest to lowest host extinction risk (as estimated by Cheung et al. 2005) 30 , parasite species richness, modularity contribution, and nestedness contribution. Many of these variables are related, which also makes the extinction orders related (see Figure S4 for rank correlations between all species removal scenarios). After each host removal, nestedness (NODF), modularity (Q), and parasite diversity (PSR; total number of parasite species remaining) were determined. Hosts were removed in random order for 1000 simulations to provide appropriate bounds of expected trajectories for the decay of nestedness, modularity, and parasite diversity. To measure the magnitude of the effect of each of the four species removal scenarios, we compared the decay curves of our targeted removal scenarios to a null expectation generated by 1000 random extinction simulations. From these simulations, we quantified the divergence of our targeted extinction scenarios to random simulations. To do this, we first plotted the decay in parasite diversity, modularity, or nestedness as a function of the number of hosts removed for each of our targeted and random simulations. Then we quantified the decay area (DA) between the curve generated from a targeted extinction scenario and the mean curve generated from random extinction simulations in terms of parasite diversity (DA PSR ), nestedness (DA C ), and modularity (DA Q ). This measure is almost identical to how Evans et al. (2013) 39 defined robustness, except they considered robustness the area under the decay curve to the origin instead of to a mean random simulated extinction scenario.
We also calculated DA values for each random simulation to the mean random simulation, thus creating a distribution of possible decay area values obtained from completely random extinction simulations. We treated this area as an estimate of the magnitude of difference between a particular host extinction scenario and random extinction simulations, which allows us to compare random extinction scenarios to our targeted extinction scenarios using z-tests. We performed all analyses in the R Statistical Computing

Results
Host traits associated with network structure. Hosts contributed differentially to network properties, resulting in many hosts having a negligible effect on nestedness or modularity, while a small number of hosts had large positive effects on network structure measures. Contrary to the idea that modularity and nestedness have opposing influences on network structure (i.e. should be negatively correlated), we found a significant positive relationship between the two measures (r s = 0.485, p-value < 0.0001).
Host-specific contributions to nestedness and modularity were most associated with parasite species richness based on the relative contribution (RC) values obtained from the boosted regression tree analysis (Fig. 1). Host family was the second-most important variable to both nestedness contributions (RC = 40.28) and modularity contributions (RC = 19.94) of host species, indicating a potential role of host phylogeny in determining the host species relative impact on network structure. The biomass of fish captured per host species had a marginal influence on host-specific contributions to nestedness. The remaining variables contributed little, each having a relative contribution value of less than 5 units (out of 100). To examine the predictive accuracy of our models, we regressed model predictions against actual values of nestedness and modularity contribution. Regression plots of model predicted versus actual values for nestedness and modularity contribution are provided in the Supplementary Materials ( Figure S3). Root mean squared error is another commonly used scale-dependent measure of predictive accuracy. Nestedness contribution (RMSE = 1.06) and modularity contribution (RMSE = 1.05) both had relatively low RMSE values, slightly over one unit difference between predicted and actual values. Both biomass and abundance are variables defined in terms of catch per unit effort, which makes both variables closely related to the dominance of a species in a community, and to sampling effort, as larger values would indicate species more frequently examined for parasites. However, sampling effort was not a strong contributor to either measure of network structure (C i and Q i ), suggesting the importance of biomass to C i could relate more to host biology and community composition than sampling effort.
Effect of nonrandom host extinctions on secondary extinctions. Sequential host removal scenarios based on extinction risk, parasite species richness, nestedness and modularity contributions resulted in significantly larger decay areas (DA) in parasite diversity(DA PSR ), nestedness (DA C ), and modularity (DA Q ) (Figs 2 and 3) compared to random simulations. The exception to this was host removal based on estimated extinction risk, which resulted in significantly larger decay area when looking at parasite diversity, but not nestedness or modularity. In other words, parasite diversity, but not nestedness or modularity, declined at a faster rate when hosts were sequentially removed based on their estimated extinction risk (Table 2).
Host removal based on parasite species richness values resulted in larger decay area values than removal of hosts based on their estimated nestedness and modularity contribution ( Table 2 and Fig. 3), suggesting both a strong link between parasite species richness and species level contributions to network structure (see Supplementary Material for rank correlation matrix) and also that nestedness and modularity contributions may change markedly as hosts are removed from the network. In other words, Scientific RepoRts | 5:13185 | DOi: 10.1038/srep13185 a policy that re-calculates nestedness and modularity contribution after each host removal may converge upon on a similar or better host removal strategy than removal based on parasite species richness.
Another interesting result is what happens to modularity values with any (random or non-random) host removal. Modularity values actually increase until a critical point, after which they go to zero. This is likely a result of the modularity calculation itself. We demonstrate this in the Supplementary Materials ( Figure S6) by numerically simulating sequential host removal in randomly constructed bipartite networks. Understanding if this increase in modularity is a statistical artifact, or an important component of complex networks is an area for further research.
There are many host removal strategies not examined here. For instance, host extinction may be related to their native or non-native status, or their abundance within the community. While these analyses may be outside of the scope of the current manuscript, we have examined both in the Supplementary Materials ( Figure S5 in particular).

Discussion
Declining host biodiversity could result in secondary extinctions of obligate parasites, especially parasites endemic to a small subset of host species. These secondary extinctions, in turn, may result in further extinctions as demonstrated recently in a host-parasitoid system 40 . Coextinction may destabilize networks 41 , resulting in faster losses in biodiversity and network structure than expected under random extinction scenarios. Here, we find that the differential impact a host has on network structure and parasite diversity can be partially explained by host traits. First, we found that hosts contributing strongly to nestedness also contributed strongly to modularity, partially reconciling two opposing network properties related to network structure 42 , though more research is necessary to determine the true relationship between these network properties. Secondly, we found that both network structure measures were related to host family affiliation and parasite community richness, suggesting that host phylogeny and the number of interactions a host has in the network play a role in determining network structure. Also, fish biomass was marginally important in predicting nestedness contributions, which may relate to biased sampling, or to host body size, as fish biomass combines aspects of individual fish mass and relative abundance. Taken together, these findings suggest that species roles in antagonistic networks may be phylogenetically conserved or influenced by host traits (e.g. body mass). Lastly, we found that random extinction simulations may underestimate the rate of decline in parasite diversity when compared to extinction based on realistic extinction risk measures, and that extinction orders based on parasite species richness, and network structure contribution measures tended to result in much faster declines in parasite diversity and nestedess relative to random extinction simulations.
Previous studies have argued that topological measures of network structure may confer stability 14,17 . For instance, nestedness 43 and modularity 17 have both been related to network persistence (proportion of surviving species after recovery from perturbation), and resilience (time until community returns to equilibrium after perturbation). Further, resilience, robustness and structural stability (all slightly different components of stability) have previously been related to nestedness 14 . These topological measures are appealing proxies for network stability as they are easy to measure and do not require time-series data to calculate. Previous studies have acknowledged that some species are more important for network stability than others 8,9,44 . Our decay area (DA) statistic can be considered a measure of network robustness 39 , thus relating topological properties of networks (nestedness, and modularity) to dynamical properties (stability).  Host removal based on parasite species richness (purple) resulted in the fastest declines of network metrics examined, followed closely by host removal based on modularity contribution (yellow) and nestedness contribution (red). Host removal by estimated extinction risk (green) did not differ from random extinction in terms of nestedness and modularity losses, but parasite diversity did decline faster relative to random simulations. While the importance of host family to species level contributions to nestedness and modularity could be viewed in an evolutionary context, it is potentially a result of the variation in parasite species richness seen in the fish taxa examined. For example, the family Anostomidae (specifically four native Leporinus species) had the highest positive values of C i and Q i and also the highest parasite species richness relative to other fish taxa. While the importance of host family may be conflated with parasite species richness, it is also possible that host family represents a coarse measure of phylogenetic relatedness, such that host family can be considered a proxy for evolutionary relatedness. Alternatively, this could simply be a variable related to unmeasured trait variation. Rank correlations between extinction scenarios show support for the relationship between species contributions to network structure and parasite species richness ( Figure S4). Overall, this suggests that parasite species richness may be the most important determinant of network structure and parasite diversity, even when compared to metrics that calculate species-level contributions to nestedness and modularity. There are a number of possible explanations for this.
First, nestedness and modularity contribution metrics were not dynamically updated as species were removed. This was intentional, as the goal of this study was to determine if a snapshot of a full network could capture node importance as the network disassembles. This result makes intuitive sense, as the removal of host species with species-rich parasite communities may be more likely to result in co-extinction of parasite species. This finding, despite its appealing simplicity may not be true for some networks, as other aspects of network structure, host biology and community dynamics may also be important for network stability 45 . Second, it may be the case that nestedness and modularity contribution measures simply cannot capture a host species true contribution to network structure, such that more simple metrics like node degree (i.e., parasite species richness) may outperform more complicated metrics. Lastly, the calculation of nestedness and modularity contributions is based on null model randomizations of each host in the network. That is, nestedness and modularity contribution measures are based on the idea that the host is still present, but the interactions of that host are randomized. Statistically, this approach makes sense, since some network structure metrics are sensitive to matrix size or fill. However, this formulation is conceptually different from a host removal, which may explain why removal strategies based on nestedness and modularity contributions did not result in the fastest declines in network structure.
Through simulated sequential host extinctions, we found that parasite diversity declined at a slightly faster rate when extinction order was based on host extinction risk than compared to random extinction order, but only after around 40% of hosts were already lost from the system. This suggests that the loss of host biodiversity based on realistic estimates of extinction risk may result in a greater loss of parasite diversity through secondary extinctions than is currently expected, but only in highly perturbed systems subject to large losses in host diversity. This result is comparable to Davies and Yessoufou (2013) 46 , in which they found that phylogenetically clustered extinctions result in disproportionate loss, compared to random extinction scenarios.
While parasite diversity and nestedness declined precipitously in our simulations, modularity actually increased slowly in our extinction scenarios relative to random extinction simulations. This finding may be a function of removing hosts with high parasite species richness as this may fragment the network into small groups of specialist parasites, which would increase modularity. In a limiting case, a host-parasite network in which each host is parasitized by a single host-specific parasite would be perfectly modular. However, modularity increased even in random extinction simulations, suggesting one of two things; either modularity statistics are biased by matrix size or host removal from a network may actually increase modularity until a critical point is reached before the network collapses entirely. We address the relationship between matrix size and modularity in the Supplementary Materials ( Figure S6).
Biodiversity loss is a pressing concern, with recent focus on the impact of host extinctions on secondary extinctions of obligate parasites 6 and the importance of parasites as indicators of ecosystem health 47,48 . A recent investigation into the shape of the relationship between free-living species richness and parasite richness under sequential host removal simulations resulted in a myriad of possible trajectories for the decay of parasite diversity 6 . While we acknowledge that host diversity loss may influence parasite diversity idiosyncratically, our simulations find a consistently greater loss in both network stability and parasite diversity relative to random simulations. Understanding the network properties that alter decay  trajectories (e.g. degree distribution) may provide a better understanding of the structure and stability of ecological networks 49 . Lastly, the integration of network structure measures into studies of network stability may yield insight into species-level characteristics that promote biodiversity and network structure.