In silico experiments

We ran digital experiments on the artificial life platform Avida version 2.14.0 (http://avida.devosoft.org/)11. The Avida platform is a model system for investigating how co-evolution (and co-extinction) works within simple, transparent, stochastic rules11,23. The most important difference between empirical food-web assembly24,25 and artificial life simulations is that food-web biologists model natural-looking systems to explain patterns in nature, whereas artificial-life simulations are alternative ‘natural’ systems that evolve in silico26,27. Artificial life simulations generate and maintain complexity starting from a few rules (just as in natural evolution). This, in turn, generates the patterns expected in ecological and evolutionary systems.

We studied host–parasite networks for several reasons. First, although Avida can model predator–prey interactions28, and can create different trophic levels by manipulating the resource setting29, these options are unexplored. Conversely, host–parasite networks in Avida have ecological and evolutionary behaviors close to those observed in nature28, and have been already used to explain how complex features evolve30 and how parasites maintain host diversity12. Furthermore, empirical food webs have limitations, such as unequal resolution among the trophic levels due to lumping some species in broad taxonomic categories31. For instance, in aquatic trophic webs, it is common to have fish identified at the species level, and all the phytoplankton aggregated into a single category32. This creates obvious problems for interpreting robustness, especially when lower trophic levels are aggregated, leading to unrealistic assumptions such as all phytoplankton go extinct at once as if they were a single species. Additionally, there is little information from IUCN on extinction risk for non-vertebrate species. Thus, even if we obtained resolved food webs, we would have not been able to investigate their robustness to future extinction scenarios. More importantly, host–parasite networks are bipartite, which makes it easy to isolate how secondary extinctions affect system robustness, whereas in food webs, hosts can suffer secondary extinctions if their resources vanish7.

We used the Python programming language33 and R34 to process Avida output, simulate disassembly, and analyse data. Supplementary Data 1 details the co-evolved host–parasite networks and the host extinction sequences in the historical scenario.

The Avida study had three distinct phases: coevolution, assessment and disassembly. In the co-evolution phase, we ran simulations to evolve 100 complex host–parasite networks (100 being enough for hypothesis testing in past studies23,30). To make our results as general as possible, we randomized several parameters (carrying capacity, parasite virulence, resource availability, injection timing and amount of ancestral parasites) of a setting already explored in host–parasite co-evolution experiments30. In each simulation, the Avida world was a bi-dimensional grid with random dimensions between 50 and 120 host units (thus between 250 and 14,400 host habitations). Mutation rates were the same as in Zaman et al.30, but we also allowed parasite virulence to mutate throughout the experiment (starting from 1, that is, a situation in which the parasites steal all the CPU to its host). To create environmental variation, we randomly selected between one and nine resources associated with the canonical nine Avida logical operations (not, nand, and, orn, or, andn, nor, xor and equ). Similarly, we randomly associated the available resources as products of the nine tasks, with an input–output ratio randomly selected between 0 and 0.5. We set a random seed for each replicate, inserted a single host ancestor and allowed for host diversification for a random period lasting between 1,000 and 5,000 steps, at which point we injected 500–1,000 individuals belonging to a single ancestral parasite species. Both the host and the parasite ancestors could only do one of the two least complex tasks in Avida (that is, the ‘NOT’ function, where 0 is returned if 1 is consumed, and vice versa)30. Parasites in Avida are similar to free-living species in terms of structure and evolutionary processes (that is, mutation type and rate). However, parasites could not survive outside a host. Thus, when a parasite reproduces, its offspring try to infect a nearby host (like a directly transmitted, single-host microparasite). The host is susceptible only if it is uninfected, and the parasite can do at least one of its tasks. Depending on their virulence, parasites can take up to all the host’s CPU cycles. During the co-evolution phase, hosts and parasites competed, interacted and co-evolved, generating complex host–parasite networks with different structural properties. We retained the first 100 simulations where at least one host and parasite species persisted to a random end point between 105 and 5 × 105 steps (resulting in assemblages with different ages).

After stopping the co-evolution phase, we assessed host vulnerability to extinction by letting host species interact but not mutate12. In this context, host species compete, going extinct one after another, depending on their relative fitness, providing an objective way to measure host vulnerability to extinction under the conditions in which they evolved in the coevolution phase (that is, historical conditions).

In both the co-evolution and assessment phases, we monitored the species every 100 steps by recording genome, spatial position, host or parasite, genetic code and closest ancestor. The Avida documentation at https://github.com/devosoft/avida/wiki gives additional details about how to setup/run Avida simulations, and processing/interpreting Avida output.

In the assessment phase (that is, after we stopped mutations), we also recorded, for each host: (i) vulnerability to extinction, measured as h/H, with h being the order a species went extinct and H being the starting host abundance; (ii) parasite richness; and (iii) average parasite host range; that is, the average number of hosts (including the target one) used by its parasites. Then, we assessed the relationships between parasite richness and average host range (standardized by, respectively, the total number of parasite species and the total number of host species per simulation), and host vulnerability to extinction (using Spearman’s rank correlation coefficient) on each individual run, and on all the runs aggregated.

In the disassembly phase, we removed hosts one by one until no hosts remained, counting parasite species richness at each step. We determined the best-case scenario for parasites by removing, in sequence, the host with the least parasites, so that hosts with many parasites went extinct late and parasite diversity declined slowly. We determined the worst-case scenario for parasites by removing, in sequence, the host with the most parasites, so that hosts with many parasites went extinct early, and parasite diversity declined rapidly. Then, to identify the parasite assemblage robustness to historical conditions, we removed, in sequence, the most vulnerable host to extinction as quantified in the assessment phase. Finally, to identify parasite assemblage robustness to a hypothetical environmental change, we randomized the host-removal sequence.

For each removal treatment, we averaged 100 replicates where we randomized ties (for example, the host extinction order associated with the same number of parasites in the best- and worst-case scenario, or having the same historical vulnerability or the same maximum complexity). In the random scenario, we simply averaged 100 random sequences. We quantified parasite assemblage robustness as the area under the host diversity versus parasite diversity curve (rescaled as proportions)7.

We considered the extent that various assumptions might affect our results. Specifically, we asked how network structure, task complexity, genotype, incomplete information and time in the co-evolutionary phase affected robustness. As detailed in the following paragraphs, these factors did not alter the qualitative findings.

Because robustness can be sensitive to network structure35, we controlled for network structure by applying each treatment to the same network. Nonetheless, our replicate networks differed in structure and that could lead to differences in robustness among them, adding to the variation we see in robustness within a treatment. To investigate whether and how network structure can affect robustness, we did pairwise comparisons between robustness, and basic network properties such as the number of nodes and edges, the connectance, nestedness/overlap and modularity36. Most networks tended toward either a nested or a modular structure or both (Supplementary Fig. 2), whereas segregation (that is, the tendency of species to minimize overlap in partners) never emerged, which might support the idea that sharing interacting partners could promote network persistence. Usually, network structure (Supplementary Table 1) did not affect robustness, explaining, at most, 10% of the variation (Supplementary Table 2).

We investigated whether task complexity (see Table 1 in Lenski et al.23) affected parasite assemblage robustness by replicating the disassembly experiments on the 100 networks and removing hosts in decreasing or increasing order of complexity. We measured host complexity as equal to its most complex task. Complexity did not correlate with robustness (P=0.21). Moreover, the robustness measured using decreasing complexity was not distinguishable from that obtained from random removal (P=0.56), whereas robustness measured using increasing complexity was only slightly different from the random removal (P=0.047; see also Supplementary Fig. 3).

To identify and disassemble networks, we classified organisms (both host and parasites) into taxonomic units (hereafter ‘species’) by their phenotype, that is, their ability to do particular tasks30. However, to assess whether evolutionary convergence affected the results, we replicated all the disassembly experiments by classifying taxonomic units by genotypes. We found consistent patterns suggesting that evolutionary convergence was not biasing our findings (Supplementary Fig. 4).

Empirical networks suffer from incomplete information that underestimates generality (for example, the more we study networks, the more complex they appear, with increasing redundancies). To better compare the simulations with the analysis on the (undersampled) empirical networks, we asked what simulations would have looked like if we had incomplete information. We replicated all the analyses by eliminating 10, 20, 30, 40, 50, 60% of the host–parasite interactions, finding that data gaps only slightly underestimated parasite assemblage robustness (Supplementary Figs 5–8).

To investigate how soon parasite assemblage robustness emerges through the co-evolutionary phase, we measured a random parasite assemblage’s robustness to historical and novel perturbations every 1,000 steps, from 10,000 to 100,000 generations. Figures 1 and 4 show that robustness to historical conditions evolved relatively early, but there is a tradeoff with robustness to novel conditions.

Tests on empirical host–parasite networks

As a companion to the Avida experiments, we evaluated our predictions using all fish–parasite records available from FishPEST (http://purl.oclc.org/fishpest)13, and 16 large host–parasite networks, built by combining all the records available from the host–parasite database of the Natural History Museum of London (http://www.nhm.ac.uk/) for different combinations of host and parasite taxa (amphibians, birds, mammals and reptiles versus acanthocephalans, cestodes, nematodes and trematodes). These networks have already been used to estimate global parasite richness14.

Both datasets were filtered by including only hosts present in IUCN redlist. After this filtering procedure, the FishPEST network included 16,681 associations between 1,696 fish species and 7,555 parasite species belonging to different taxa (acanthocephalans, cestodes, monogeneans, nematodes, trematodes), whereas the Natural History Museum of London data included 29,275 associations between 2,751 host species and 9,638 parasites species.

We disassembled these networks, by taxon, by removing hosts as in our experiments on digital networks. For all host groups, we could estimate future vulnerability from IUCN risk categories for some species. Thus, for both the vertebrate and the fish–parasite networks, we simulated 100 worst-case, best-case and novel scenarios (randomizing the whole order of extinctions in the random scenario, and only ties in the others). For the fish–parasite network only, we simulated 100 disassembles under historical conditions using fish intrinsic vulnerability to extinction (www.fishbase.org)15 as a proxy for historical vulnerability. Such measures, which take into account fish life history and ecological characteristics, have already been used as a proxy for fish intrinsic extinction risk in studies dealing with co-extinctions8,18.

Sensitivity analyses

We considered the extent that various assumptions might affect our empirical results. For instance, we calculated robustness as if all parasite species had simple life cycles. However, many parasites require two or more host species in sequence, and this reduces robustness7. Another aspect that reduces robustness in food webs (absent from our disassembly) is that hosts can suffer secondary extinctions if their resources vanish7. However, because the overestimation applies to all the scenarios, we expect no bias in the results. The following paragraphs describe how we determined that reducing the FishPest data filtering with IUCN records, ties in the IUCN vulnerability rankings, incomplete parasite information and sampling-effort bias did not alter the qualitative findings.

Although filtering the FishPEST data set by IUCN records reduced the data set, our analysis on Avida partial networks suggested that incomplete information does not introduce biases, instead underestimating parasite assemblage robustness. We found the same pattern by replicating the fish host removal according to intrinsic host vulnerability on the complete FishPEST network (including 33,426 unique interactions between 12,762 parasite species and 4,091 host species). Supplementary Figure 9 plots parasite diversity versus host diversity for, respectively, the complete FishPEST network, and that filtered by IUCN data (that is, the same curve shown in Fig. 5).

A potential problem with comparing the empirical data is that the IUCN categories lead to more ties in the vulnerability ranking than occurs for the continuous FishBase vulnerability measure15. To investigate whether ties affected the results, we replicated the fish host disassembly by grouping fish into five intrinsic vulnerability categories (corresponding to the five intervals in the FishBase vulnerability scale). The resulting historical disassembly curve was almost identical to the one obtained using continuous vulnerability values (R2=0.99; regression line: slope=1.0; intercept=−0.000004).

Furthermore, the species that are most under-sampled for parasites are either rare, or limited in range. According to our main hypotheses and findings, rare species are likely to be used by few generalist parasites. Thus, their addition to the data sets would not have a strong effect on parasite persistence. Similarly, missing some host records for widespread generalist parasites would affect the overall robustness little, because these species are already likely to persist. Even if the estimate error was appreciable, we do not expect it to affect the comparison among treatments. The biggest concern about bias for the empirical data would be if increasing host sampling effort disproportionately sampled vulnerable hosts and found more specialist than generalist parasite species. Such a pattern could make it appear that invulnerable hosts had more specialist parasites when they were simply sampled more. Fortunately, sampling effort (measured as published parasite studies per fish species) did not correlate with vulnerability (r s =−0.015) or specialist to generalist ratio (r s =−0.049).

Data availability

The 100 co-evolved host–parasite networks and the host extinction sequences in the historical scenario are available in Supplementary Information. Host–parasite records used in the analyses are available from FishPest (http://purl.oclc.org/fishpest) and from the London Natural History Museum host–parasite database (http://www.nhm.ac.uk). Fish vulnerability scores can be obtained from FishBase (http://www.fishbase.org), while conservation status information for both fish and terrestrial vertebrates can be retrieved from the IUCN website (http://www.iucnredlist.org).