Multilayer networks reveal the spatial structure of seed-dispersal interactions across the Great Rift landscapes

Species interaction networks are traditionally explored as discrete entities with well-defined spatial borders, an oversimplification likely impairing their applicability. Using a multilayer network approach, explicitly accounting for inter-habitat connectivity, we investigate the spatial structure of seed–dispersal networks across the Gorongosa National Park, Mozambique. We show that the overall seed–dispersal network is composed by spatially explicit communities of dispersers spanning across habitats, functionally linking the landscape mosaic. Inter-habitat connectivity determines spatial structure, which cannot be accurately described with standard monolayer approaches either splitting or merging habitats. Multilayer modularity cannot be predicted by null models randomizing either interactions within each habitat or those linking habitats; however, as habitat connectivity increases, random processes become more important for overall structure. The importance of dispersers for the overall network structure is captured by multilayer versatility but not by standard metrics. Highly versatile species disperse many plant species across multiple habitats, being critical to landscape functional cohesion.


Supplementary Information
Supplementary Methods DNA barcoding of seeds DNA was extracted from seeds using the DNeasy Plant Mini Kit (Qiagen) following the manufacturer instructions. The DNA samples were subjected to an additional step of purification with Phenol: Chloroform: Isoamyl alcohol. DNA was re-suspended in 40 μl of elution buffer and kept at 20 ºC. Two chloropastid loci (the psbA-trnH intergenic spacer, and the trnL intron and trnL-F intergenic spacer) were amplified using a Hot Start Taq Master Mix (QIAGEN) as described in 1 . Amplification was performed in 25 μl containing 1 μl of DNA and 1 μl of each primer. Conditions of the PCR were as follows: 95ºC (15min); 94ºC (1min); then 30 cycles for trnL-F and 35 cycles for psbA at 94ºC (1min)/ 50ºC (1min)/ 72ºC (1min), and a final extension at 72°C (10 min). The PCR products were purified using ExoSAP-IT (Affymetrix), and sequenced in a Sanger ABI 3730xl at GATC Biotech (Germany). The sequences were compared with the available online databases using BLAST 2 . The species were identified based on the best BLAST matches and the list of plant species known for the Gorongosa National Park.

Multilayer modularity
Modularity is a structural pattern of interactions between nodes of a network whereby a group of speciesa module or community, interact more frequently than expected among them than with other groups of species 6,7 . A multilayer approach to modularity allows the identification communities that span across multiple layers of the network, which can be important to the structural unity of the whole network 8 . We used a modularity quality function that uses a "generalized Louvain" method to community finding 9,10 . The Louvain method for the identification of communities progresses in two iterative phases: in the first phase, all nodes are considered one-by-one and assigned to a specific set of nodescommunity, until a configuration is reached that maximizes the modularity quality function. In the second phase, the communities previously found are now used as nodes of a reduced network, and the same procedure is repeated until no further increase in modularity is detected 11 . This is a popular locally-greedy method for modularity-optimization as it is fast and delivers reliable results 12,13 .
where A ijs is the weight of the intra-layer edge between nodes i and j within layer s; C jsr is a tensor element giving the weight of the inter-layer between node j and its replica on layers r and s (given the categorical nature of the multilayer coupling in spatial multilayer networks all values C jsr > 0, and it is assumed to be equal for any inter-layer coupling, C jsr = ω); γ s is the resolution parameter for layer s; k is and d js are the degrees of plant i and dispersers j within layer s, respectively; m s is the total edge weight of layer s; g is and g jr are the set of nodes forming the communities that contain the nodes-layer (i,s) and (j,s), respectively; the Kronecker delta between indices x and y is denoted as δ xy (this will be 1 for x = y and 0 for x ≠ y), and 2µ = ∑ ijs A ijs 9 .
The "generalized Louvain" methods requires the specification of two parameters: the resolution limit γ; and the inter-layer coupling ω. The resolution limit γ defines the detail to which the network will be resolved into communities, and can be seen as the importance given to the null model relative to the empirical network 12 . We used the default resolution parameter value of γ = 1 9,12 . The choice of the coupling parameter ω is a matter of intense investigation, and takes a value of either 0 or ω 12 . When ω = 0 it is equivalent to optimizing the modularity for each layer independently, where any node never belongs to the same community across the different layers, i.e. communities are not persistent across the multilayer network. If however ω > 0, and as it increases, nodes are less likely to belong to different communities, which tend to span across the different layers of the network, and can assume different values of each pair of layers depending on the importance of the coupling between those pairs of layers 9,12 .

Versatility
To assess the importance of nodes to the structure we calculated centrality for each node accounting for the multilayer nature of our network, defined by the animal-plant interaction in each of the habitats of Gorongosa. This allows to identify the most important nodesversatile species, in our system 15 . We used a widely used measure of centrality based on Google's PageRank 16 , which is a random walk centrality measure corresponding to the path taken by a walker moving between adjacent nodes, with the importance of each node being calculated recursively by the sum of the importance of all nodes connected to it. PageRank centrality was extended to the case of multilayer networks by allowing "teleportation" of nodes between any layers of the network 15 .

Multistrength
Node multistrength measures the strength of a node as the combined weight of its connections, across the different layers of a network 17,18 , and expresses the importance of a node to the community of nodes with which it interacts in the multilayer network.
Two concepts are important to understand multistrength, namely: multidegree and multilink. Multidegree is the number of links in which a node participates, and it is an extension of node degree for monolayers 17,18 . A multilink is defined as the set of links