Composition, structure and robustness of Lichen guilds

Symbiosis is a major engine of evolutionary innovation underlying many extant complex organisms. Lichens are a paradigmatic example that offers a unique perspective on the role of symbiosis in ecological success and evolutionary diversification. Lichen studies have produced a wealth of information regarding the importance of symbiosis, but they frequently focus on a few species, limiting our understanding of large-scale phenomena such as guilds. Guilds are groupings of lichens that assist each other’s proliferation and are intimately linked by a shared set of photobionts, constituting an extensive network of relationships. To characterize the network of lichen symbionts, we used a large data set (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=206$$\end{document}n=206 publications) of natural photobiont-mycobiont associations. The entire lichen network was found to be modular, but this organization does not directly match taxonomic information in the data set, prompting a reconsideration of lichen guild structure and composition. The multiscale nature of this network reveals that the major lichen guilds are better represented as clusters with several substructures rather than as monolithic communities. Heterogeneous guild structure fosters robustness, with keystone species functioning as bridges between guilds and whose extinction would endanger global stability.

Lichens are symbiotic organisms composed of a fungus (mycobiont), one or more photosynthetic partners (photobionts, typically algae or cyanobacteria see Fig. 1a, b) and other microbial species 1,2 . In the wider ecological context, lichens provide several services that are essential for ecosystem functioning: from weathering of rocks increasing the bioavailability of minerals 3 to carbon and nitrogen fixation 4 . By virtue of the wildly different metabolisms of photobionts and mycobionts, lichenization provides new biological traits that can enable both partners to colonize a wider range of environments 5,6 , including extreme 3,7 , polluted 8,9 or anthropogenic ecosystems 10 . In consequence, lichens are more than the sum of their constituent symbionts, emphasizing the relevance of non-fraternal organismality as a source of evolutionary innovation 11,12 . Symbiosis is a natural hot spot for diversity and innovation 13,14 , the evolutionary potential of a lichen symbiont must take into account the extraordinary repertoire of traits gained by acquiring a new partner. However, it is widely acknowledged that specificity of the mycobiont-photobiont association as well as the spatial distribution of their component species largely drives the formation of new lichen organisms 6,[15][16][17][18] . Indeed, many symbionts are quite stringent in the partnerships they form, which are often determined by biophysical gradients [19][20][21][22][23][24] or substrate preferences 18 . Selective mycobionts do not always use the entire niche of their photobionts and are therefore sometimes restricted to certain climatic conditions. It has also been suggested that highly specific mycobionts are typically restricted to few environments and display limited ecological range 25,26 . Conversely, mycobionts capable of colonizing different habitats often are less stringent in their partnerships, interacting with various available photobionts 13,17 .
Collectively, the set of natural associations between photobionts and mycobionts defines a network of symbionts 27,28 , where each interaction corresponds to a singular lichen species (Fig. 1c). It has been proposed that within this network lichen species organize in communities known as photobiont-mediated guilds 29 . Guilds are groups of lichen species that are ecologically connected by sharing one or more photobionts 29,30 (Fig. 1d). Fungal species in a guild can benefit each other by propagating a common set of partners, driving the establishment of connected lichen species into new or marginal habitats 31 , while competing for space and resources 32 .
Much of our current understanding on photobiont-mycobiont partnerships comes from studies involving a few lichen species or geographical areas 17,19,23,24,[33][34][35][36][37][38] . However, mycobiont-photobiont partnerships do not happen in isolation, they are part of a larger web of interactions supporting the assembly and maintenance of communities, and understanding them requires a systemic approach. Many species interaction networks are classified as nested or modular. A modular network is made up of sparsely linked clusters of dense subgraphs [39][40][41] . These communities or modules may arise due to evolutionary and environmental constraints 42 , or they may represent

Results
Defining the global photobiont-mycobiont association network (PMAN). Previous research has focused on lichenization among few species of mycobionts and photobionts, providing an insufficient understanding of symbiotic interactions at the largest scales. For example, symbiotic ties could be have been shaped by the presence of additional interactions in a larger community. In this context, network techniques have been particularly successful in the study of mutualistic 39 and antagonistic 46 networks from small to large spatio-temporal scales 48,49 . Here, we recreate the largest network of lichen symbionts to date using the full data set of symbiont pairings assembled in a recent meta-study by Sanders 52,53 . Figure 3a shows the distribution of modularity values for an ensemble of bootstrap randomizations (histogram with shaded region) compared to the real data set average modularity (dashed vertical red line). This suggests that the lichen network is highly modular, more than the expected value for the null models ( p < 10 −5 ), consistently with a significant decrease in nestedness ( p < 10 −5 , see Fig. 3b and Table 2). Following standard analyses of data completeness 54,55 , we studied the modularity and nestedness of subsamples using half the dataset, i.e. real networks and their randomizations with half the number of observed network links (see Fig. 4). We discovered that modularity and nestedness are maintained in both circumstances, implying that the presented patterns are robust and independent of sampling depth.

Figure 2.
Photobiont-mycobiont association network. Bipartite network representation of the full dataset analysed in our study. Interactions in this network involve two different types of nodes: blue nodes correspond to mycobionts, and yellow nodes indicate photobiont species. The network consists of 34 isolated components, the largest of which has average degree �k� = 2.524 (see Table 1). Network layout was automatically generated with the FMMM algorithm (see "Materials and methods"). been traditionally linked to photobiont identity (or "photobiont-mediated guild"), and more precisely, at the genus level 29,33,54 . Following this approach, Fig. 5a shows the PMAN where each mycobiont is labelled using the genus of their closest photobiont (see "Materials and methods"). The labelled network comprises several contiguous clusters containing related mycobiont species, i.e., belonging to the same genus. Here, the network appears to display a higher density of connections within each of the guilds than with outside nodes, which suggests that guilds display a common structural signature based on modularity.
We compare topological modules to well-known guilds to study the link between taxonomically-specified guilds and modularity (see "Materials and methods"). In the PMAN there are 56 taxonomy-defined guilds and 140 topological modules. We find that topological modules are statistically smaller than the taxonomy-defined guilds, with an average size of 7.61 species per module versus 26.11 species per guild ( p < 10 −4 , see Fig. 5b). Smaller guilds have a strong matching with topology-predicted modules, whereas larger guilds include several smaller modules embedded within them (see the highlighted cases of Asterochloris and Trebouxia). This is computed using Jaccard's similarity indices 56 between each guild (defined by taxonomy) and each module (defined by topology). A guild or a module here are the set of fungal and photobiont species that have a given label. There   www.nature.com/scientificreports/ are 21 guilds with average Jaccard similarity less than 0.5, indicating a significant discrepancy between topology and taxonomy, while the remaining 35 have very similar communities (Fig. 5c). The lack of matching between taxonomic and topological classification implies the existence of other structural patterns beyond the mesoscale defined by guilds. Empirical support to this hypothesis is found in the internal structural diversity of two major guilds: Coccomyxa and Trebouxia. These guilds are well established in the literature 29,32 and provide a reliable microcosm whose individual components been extensively characterized, making them an ideal test-bed for our approach. Figure 6a shows a schematic representation of guilds proposed by Rikkinen 29 , while Fig. 6b, c, on the other hand, shows the sub-networks of the same guilds with the most current information. These representations intuitively show distinct properties, but the network approach allows us to quantify these differences more precisely.
Guild interconnectivity matches major photobiont groups. At a bigger scale, it is important to assess whether there the interconnectivity among guilds is structured or not. Are there any non-trivial patterns beyond the level of guild organization? We build a guild interaction network (GIN) where nodes represent the major guilds (consisting of at least 5 species) and edges connect each pair of guilds that share at least one mycobiont species (see "Materials and methods"). Figure 7 depicts the adjacency matrix (a) and the circular layout (b) for the GIN. Using adjacency as a surrogate for distance in this network, we can estimate node similarity and lay out guilds using the fastcluster algorithm 57 . The photobiont phylogeny can be partially recovered in the GIN: several large green algae guilds cluster together (from Heveochlorella to Asterochloris) followed by smaller modules containing some green algae, Xanthophytes and cyanobacteria located in its own domain. In Figure 7b, the network representation allows us to show additional information: edge width is proportional to the number of shared species amongst guilds, while node size is proportional to the number of species in each guild. Guilds are colored according to photobiont type, green for green algae, blue for cyanobacteria and yellow for yellowgreen algae. In agreement with our prior clustering analysis, every cyanobacterial guild groups together except for Nostoc. This implies that the large-scale organization of guilds is not random, but rather the result of biological traits or evolutionary constraints that bring physiologically similar species closer together in the network.
Guild and species contribution to network robustness. We investigate the PMAN's ability to withstand various types of perturbation while also providing some topologically-based insights to preservation efforts through the identification of keystone species and guilds. The degree of network fragmentation caused by species extinction is quantified using global efficiency, which reflects how costly it is to convey information across nodes. Global efficiency falls as the distance to a given node increases, eventually approaching zero in a fully disconnected system (see "Materials and methods"). Figure 8a shows the effects of different strategies of species removal on the global efficiency of the PMAN: random removal of species (blue), removal based on degree (green) and removal based on centrality (yellow). The PMAN is particularly resistant to random component failure; even after removing 120 nodes at random (about 10% of the network), the overall efficiency of the system is barely affected. This is consistent with many observations in heterogeneous biological networks 58   www.nature.com/scientificreports/ but most species have few connections and their extinction has a much more limited influence on the system's structure. In contrast, the network is particularly vulnerable to targeted node removal, whether based on degree or node centrality. While topological network properties are unlikely to inform ecological damage, this sort of study might help guide conservation strategies that seek to prevent fragmentation and the buildup of ecological damage. Figure 8b shows a simulated cascading extinction event, in which random species disappear and with a given probability this effect is propagated to neighboring species. When the chance of propagation increases, the amount of affected species increases, reaching the totality of the system when the probability is 1. This type of analysis is carried as a comparison between the the real data set (blue) versus an ensemble of bootstrapped networks (orange) acting as negative controls. On average, the real data set is more robust to cascading extinctions   Figure 7. Global connectivity between photobiont genera. Using the photobiont genus as the foundation for lichen guilds, we reconstruct the guild network's adjacency matrix (a), where white squares represent at least one mycobiont shared by photobiont genera. Photobionts are sorted using their pairwise distances and fastcluster algorithm (see "Materials and methods"). In (b) we show the corresponding network representation, where link widths are proportional to the amount of shared mycobionts between photobiont genera, node sizes are proportional to the number of species (both photobionts and mycobionts) belonging to that guild and and are colored according to their group: green for green algae (Chlorophyta and Charophyta), blue for cyanobacteria and yellow for yellow-green algae (Xanthophyta). Network visualizations were generated with OGDF library (see "Materials and methods"). Some guilds and species have a greater influence than others in providing global resilience due to their placement within the network as well as their local characteristics (Fig. 8c). Of particular interest are the guilds Scytonema and Gloeocapsa, which act as bridges to the major cyanobacterial guild Rhizonema. In addition, species within the major guilds Elliptochloris, Symbiochloris and Chloroidium have relatively high centralities. Remarkably, these are fragmented guilds (they do not form continuous clusters in Fig. 5a) and their constituent species often lie at the intersection between other major guilds in the network.

Discussion
At the core of the lichen symbiosis is the association between two unequal partners, a fungus and a phototroph. The flexibility of this association is key to understand the resilience of lichen species to polluted and antropogenic ecosystems 10 and their capacity to adapt to new environments 6 . Understanding and documenting single partnerships has been a major focus of lichen research. However, mycobiont-photobiont partnerships do not happen in isolation, they are part of a larger web of interactions that shapes their chances of formation and success. While much effort has been invested in characterizing individual associations, the global network of interactions that maintain the lichen symbiosis remains unknown. Using the most comprehensive data set of mycobiont-photobiont partnerships to date, we addressed the multi-scale nature of lichen symbiosis from a network perspective.
Ecological networks have emerged as a powerful tool to formalize biotic interactions, including mutualistic, antagonistic and more complex relations 59 . This characteristic makes them particularly useful to study the so called photobiont-mediated guilds 29 : communities of lichens that share a common set of photobionts and facilitate each other's propagation, but also compete for space and resources. Lichens depend on a set of coupled ecological interactions involving species dispersal, facilitation and competition, which are in turn shaped by guild structure 32 . Mycobionts within a guild are able to promote the establishment of one another into new habitats by extending the ecological range of their partners 33 . For example, spore-producing lichens that require a compatible photobiont upon reproduction are indirectly facilitated by the asexual lichens within their guild already established in that environment 30 .
Previously, guilds have been defined by photobiont taxonomic identity, but this characterisation has been constrained by the scope of the data sets, which have often included a limited amount of lichen species or geographical locations 17,23,[32][33][34][35][36] . An aggregated network perspective allows us to tackle the structural signature of lichen guilds and provides a solid foundation to address how the general patterns of association in the lichen symbiont network affect its robustness.
Our analysis shows that the PMAN displays a statistically significant modular and non-nested organisation when compared to the null model. This is consistent with previous observations of Peltigera lichens 38 . Furthermore, we compare the processes of guild allocation based on photobiont taxa and topologically defined guilds, testing pre-established assumptions regarding guild structure and composition. Beyond the proposed guild mesoscale, we can identify other relevant scales in lichen symbiosis. Small taxonomy-defined guilds usually match with topological modules, and are typically composed by a single photobiont linked to a few mycobionts. Instead, larger guilds like Asterochloris and Trebouxia are rather complex communities with many interacting sub-modules. This newly reported scale in the PMAN suggests that additional underlying ecological or generative constraints (such as habitat range or trait-dependent associations) may play an important role in shaping guild structure.  www.nature.com/scientificreports/ According to the GIN analysis, trait-dependent associations may also inform large-scale interactions between photobiont genera. We found that guilds are connected in a non-trivial manner consistent with the evolutionary history of photobionts, separating algae and cyanobacteria into two identifiable clusters. A possible explanation to understand these findings might be the unique way by which some cyanobacteria are lichenized. Lichens harboring cyanobacterial and algal partners segregate them in the same organism using specific organs called cephalodia 60 . The use of specific and non-universal morphological structures for cyanobacterial symbionts points at a distinct proximal evolutionary origin for these innovations and highlights the importance of morphological novelty as a driver of the patterns observed in the GIN. Putting these two structural patterns together, the PMAN recapitulates evolutionary history at the largest scales, but is inconsistent with taxon-based guild definitions at the smallest scales.
The network approach has the potential to give insights into the underlying symbiotic associations. The nature of the lichen symbiosis, whether mutualistic, antagonistic, or somewhere in between, is a subject of debate within the community [61][62][63][64] . Some authors regard lichen symbiosis as mutually beneficial among partners 36,64 , while others propose a more complex relationship analogous to photobiont domestication by mycobionts 63,65 . Under the lens of facultative mutualism, symbiosis can be malleable, defined by additional factors like species density 66 or the quality of the external environment 62 . In harsh conditions, photobionts in association with fungi are better protected from external challenges like predation, UV radiation and extreme temperatures 67 . Yet in optimal environmental conditions there are opportunity costs to living in association; the photobiont might be better off free-living instead of transferring metabolites to its host in exchange for a protection that is not required.
Theoretical studies have used different arguments to explain the presence (or absence) of modularity and nestedness in ecological networks, including ecological stability, spatial constraints, or the strength of coevolution [68][69][70] . Evolutionary network models have shown that the topological features of an ecological network determine its resilience, stabilizing structural patterns based on the types of interactions it shows. For instance, nestedness is often the key pattern found in mutualistic networks 39,68 , while modular structures are commonly associated with antagonistic networks 42,68,69 . Modularity and nestedness are negatively correlated but they can also co-occur in sparsely connected random networks 43 . This observation highlights the remarkable structure of the PMAN, which even at low connectance values displays modular, non-nested structures. An important aspect is that nestedness can also be produced for-free by the evolutionary generative rules 71 , and it is still unclear if nestedness correlates with stability according to more recent metrics 72 . Nestedness and modularity are not exclusive to particular interaction types 73 , and they can also coexist, occupying different scales in the ecological network (e.g. modules that are themselves nested 46,74,75 ) or different dimensions of hypergraph ecologies 76 . Our results suggest the presence of non-mutualistic interactions among lichen symbionts, whether antagonism or facultative mutualism 42,68,69 .
Network structure can have a significant impact on the resilience and persistence of species and ecosystems 58,77 . How perturbations propagate in ecological networks has been the subject of numerous studies, bridging the gap between topology and population dynamics. Fragmentation of the ecological network can have a detrimental influence on effective diversity and disrupt the flows of matter, energy, and ecosystem services 78,79 . In turn, diversity has a buffering effect in accumulation of damage and the effects of perturbations in ecological systems 80 . However, not all species have the same contribution to network resilience. Removal of central or highly connected species can pose larger threats to network integrity (see Figure 8c). A modular structure also impacts the spread of cascading extinctions, confining perturbations to small groups of species 69,81 . The symbiotic network presented here benefits from both of these patterns: hub species keep the network connected even when individual species are randomly removed, and the highly modular structure increases resilience by functioning as a firewall compartmentalizing damage.
In conclusion, we found evidence of the modular signature of guilds in the dataset of aggregated photobiontmycobiont associations. However, in 25 of the 56 guilds there is lack of agreement between taxonomically defined guilds and topological modules, revealing a wealth of structural patterns operating at multiple scales. This stands in stark contrast to the largest patterns of interaction in the network, which closely follow major evolutionary groups. Our analysis shows that some guilds have a greater influence than others in terms of network robustness. Future work will have to address the biological reasons explaining the origin of the reported structural patterns in symbiont networks. Incorporation of topological patterns in theoretical models will improve our knowledge about symbiotic network robustness, which is essential to anticipate biodiversity losses and extinction cascades.

Materials and methods
Data set. To build the bipartite network of photobiont-mycobiont associations, we use the data from a recent meta-study 47 . This work compiles over 200 publications in the field of lichen symbiosis, describing natural interactions between symbionts. The data was supplemented with additional well-known interactions (Lobaria's cyanobiont partners) that were missing in the original data set. The meta-study by Sanders and Matsumoto focuses on publications after 1988 Tschermak-Woess' influential review and especially those publications providing molecular evidence of species identity. This dataset is, however, biased towards geographical locations in the western world (namely, Europe and North America). Most symbiont pairs are genetically validated through sequencing efforts and correspond to a recent wave of publications in the field, but the earliest publications in the data set (e.g. "A monography on algal culture" by Chodat, 1913) do not provide the same degree of validation. Thus, the data set contains some degenerate information in which, symbionts are identified only up to the genus level instead of the species level ( ≈ 30 % of the data). Similar patterns as those reported here were found in the analyses pertaining the non-degenerate network, as well as those limited to the largest subweb.
From the list of naturally observed associations, we can define a photobiont-mycobiont association network (PMAN) G = (P, M, E) in terms of the two disjoint sets of photobiont P and mycobiont M species and another Scientific Reports | (2023) 13:3295 | https://doi.org/10.1038/s41598-023-30357-w www.nature.com/scientificreports/ set of edges (i, j) ∈ E capturing the symbiotic association between partnering species i ∈ P and j ∈ M . The graph G can be also written as a biadjacency matrix B ij with dimensions corresponding to the size of each set of species. Any given entry in the matrix equals zero if there is no interaction between species i and j and equals 1 otherwise. From the PMAN, we can further define a guild interaction network (GIN), Ŵ = (V , I) composed of a set of guilds u, v ∈ V and their interactions I. Links (u, v) ∈ I in this network indicate whether a pair of photobionts belonging to different guilds share one or more mycobionts. We can compute the (weighted) adjacency matrix A = [A uv ] of the GIN as follows: where the membership function g i indicates what guild the photobionts species i belongs to, and δ(a, b) = 1 when a = b or 0, otherwise. The V × V square matrix A represents an unipartite network whose entries measure how many mycobionts are shared between pairs of photobiont genera.
Network visualization. Network visualizations were generated with Python 3.9 scripts, using NetworkX version 2.8 82  Taxonomical definition of guilds using label propagation. In order to establish taxonomical guilds we use a variation on the classical label propagation algorithm by Raghavan, Albert and Kumara 85 . In our bipartite label propagation algorithm, photobiont nodes are labeled with their genus, while mycobionts are left unlabeled. Mycobiont nodes then acquire their surrounding photobiont tags sequentially but in a random sequence, and each community grows isotropically. In the case that a mycobiont has neighbors in different communities, a dual community allocation is proposed for that given node. At the end of this iterative process, nodes tagged with the same labels are grouped together as guilds. This is a simple and effective method for discovering any underlying structure in the PMAN network that is led by taxonomic information.
Network metrics. Modularity. A number of metrics have been used to assess the modular structure of ecological systems, but a popular method examines the degree of modularity at the topological level 86 . Given a partition of the set of nodes given by the function c(i) indicating the module label to which the node i belongs, modularity Q calculates the difference between the proportion of links inside each module and the predicted fraction when connections are randomly rewired: where A = [A ij ] is the adjacency matrix of the (unipartite) network, 2m = ij A ij is twice the number of edges, k i = j A ij is the degree of node i, and δ(a, b) = 1 when a = b or 0, otherwise. The previous definition must be modified for a bipartite network to add the requirement that no connections exist between nodes of the same type 87 . In this case, the bipartite modularity ( Q B ) is defined as: where B = [B ij ] is the bipartite adjacency matrix, and q i = j B ij and r j = i B ij is the degree of photobionts and mycobionts species, respectively. Here, we measured the bipartite modularity of the PMAN using a custom implementation of the biLOUVAIN algorithm 88 in Python 3.9 (see link to open repository in data availability), which heuristically finds a modularity maximizing arrangement of photobionts and mycobionts into a nonpredetermined number of modules. Initially, each node is allocated to its own community and then, iteratively, the algorithm finds and applies a community allocation change (swapping the module of a single node) that yields an increase in the overall modularity. This process stops when the algorithm cannot find a community swap that further increases the global modularity beyond a predetermined threshold ρ (in our case ρ = 10 −6 ). The bipartite modularity Q B assesses how often a particular annotation of nodes into modules corresponds to interactions that are mostly inside each module ( Q B = 1 ) versus mostly outside of each module ( Q B = −1).
Nestedness. Nestedness is a global measure of the propensity of low-degree species to interact with a subset of highly interconnected species. This pattern is clearly identified by the network architecture; nestedness is a systematic arrangement of non-zero entries in the adjacency matrix. We compute the nestedness using the overlap and declining fill (NODF 89   www.nature.com/scientificreports/ value implies that certain species' interactions are a subset of other generalist species' interactions (and so the network demonstrates nesting), whereas a low value indicates clustering (which is consistent with high modularity). After determining the empirical nestedness values, they are compared against null model predictions to assess their statistical significance. Nestedness is sensitive to network size N = |P| + |M| as well as network connectivity, according to theoretical and empirical investigations. We must compare empirical data with a bootstrap technique that keeps the size and fill of the PMAN's adjacency matrix since we do not know the underlying null distribution of the test statistics.
Global efficiency. Latora and Marchiori 90 introduced a network metric of global integration that is inversely proportional to the average path length. The mean global efficiency (E(G)) of a network G can be defined as follows: where d ij is the length of the shortest path connecting any pair of nodes i and j, and N is the network size. Distances between pairs of nodes in disjointed subweb are infinite, as there is no path that can reach from the original node to the destination. An advantage of using the inverse of path length (instead of the alternative average distance) is that it allows computation of a finite efficiency value for graphs with fragmented networks. Efficiency calculations were carried out using the NetworkX version 2.8 of the algorithm for Python 82 .
Betweenness centrality. Centrality measures propose a method to evaluate node or link importance given a certain network topology. Betweenness centrality (BC) of a given node is the fraction of shortest paths that pass through that node, considered over all other possible pairs of nodes within the network. A low betweenness centrality ( BC = 0 ) means that the node is not in the path to reach other nodes, and thus has no effect in the flow of information within the network. High Betweenness centrality ( BC = 1 ) means that the node is always involved in the traffic of information in the network and its removal could have major impacts in the network's information flow, whether through rerouting of shortest paths (making them longer and thus a less well connected network) or simply fragmentation impeding further communication between sets of disjointed nodes. Here, a bipartite or unipartite implementation is essentially the same, but a consideration must be given to the sampling mechanism of node pairs. The process can be substantially sped up at the cost of reliability by reducing the fraction of pairs considered in the analysis. All calculations shown in this publication correspond to the full set of nodes, computed following the standard algorithm by Borgatti and Halgin 91 , as implemented in NetworkX version 2.8 82 .
Robustness to species loss. We tested the effects of species removal on the global efficiency using three different strategies for choosing nodes: (1) Random, (2) Degree-based, and (3) Centrality-based. For the two latter strategies node choice was carried out using a propensity vector incorporating the normalized degree sequence and the normalized betweenness centrality (calculated as described in the prior section). We also analyzed the effects of species extinction propagation in the symbiont network through a simulation of extinction cascades 92,93 . In these simulations, a node is randomly chosen from the network (both in the empirical data set and the bootstrapped data set) and an "infection" process is carried out. While a node is infected (or becoming extinct) it can affect neighboring nodes with a given probability (a parameter fixed for each simulation). If new nodes become extinct then the process continues until no new extinction events are produced. The reported results are the fraction of surviving species at the end of the extinction cascade. These simulations are carried out in 100 independent replicates (meaning a replicate for 100 different bootstraps of the original data set for the case of the negative control).

Statistical validation. Statistical significance of modularity and nestedness under sampling biases. Follow-
ing standard practices in studies of symbiotic networks 55 , we further validated the patterns reported here by performing similar analyses on randomized subsets including only half the documented interactions for data completeness 94 . We compared a set of 50 subsampled networks with the same amount of bootstrap edge randomization of these subsampled networks. The subsampled networks containing 50% of the total interactions were more modular and less nested that their bootstrap counterparts, reinforcing the perspective that the patterns reported here are meaningful and unlikely to be caused by limited sampling of the underlying interactions. Significance of differences in these metrics was quantified with a scipy t-test statistic, yielding the significance p-values of p < 10 −5 for both modularity and nestedness.
Null models of bipartite networks. Validating the statistical significance of structural patterns is a common approach in network studies (i.e., two classic examples are motif analysis 95 and community detection 96 ). The goal here is to determine when empirical patterns depart from a baseline null model that provides the expected behaviour of a system. In this paper, we employ a bootstrap randomization of the bipartite network which preserves the original degree distribution 41 . If each link connects two nodes, one of each class (photobionts and mycobionts), we break correlations in the system by untying connected nodes and connecting them at random. More specifically, we create the vectors − → M , − → P of length L (where L is the number of interactions in the PMAN) and each pair (P i , M i ) is composed of a photobiont identifier and a mycobiont identifier capable of forming a lichen species. Then, we randomize the order of the elements inside vectors − → M , − → P separately. Each