DOI: 10.5061/dryad.sj3tx96dp
title: Interactions outside local patches contribute to the compound topology of plant-pollinator networks in fragmented dune slacks

description: Plants and pollinators engage in ephemeral interactions that are essential
 for both plant reproduction and pollinator fitness. Most plant species
 interact with multiple pollinator species and, similarly, most pollinators
 visit multiple plant species, leading to complex networks of interactions.
 A better understanding of the factors that affect the structure of these
 networks is crucial for effective biodiversity conservation. While
 previous research has shown significant negative impacts of habitat
 fragmentation on plant and pollinator diversity, far less is known about
 how habitat fragmentation contributes to the overall network structure. In
 this study, we used DNA metabarcoding to construct pollen transport
 networks across nineteen dune slacks along a fragmentation gradient at the
 Belgian coast. For each network, we investigated the impact of patch area,
 connectivity, and floral diversity on network size, connectance,
 nestedness and modularity. We found a total of 3222 unique pairwise
 interactions involving 139 insect species and 261 plant species, with a
 substantial proportion of pollen coming from plants growing outside dune
 slacks, indicating frequent foraging beyond local patches. Networks were
 not significantly nested, while most were significantly modular and showed
 a hierarchical compound structure. Importantly, no significant
 relationships were found between network measures and patch area or
 connectivity, indicating that local patch characteristics had little
 impact on network structure. We conclude that plant-pollinator networks in
 fragmented landscapes with a permeable landscape matrix are able to
 maintain network characteristics, in part through interactions with plants
 outside the local patch, and that pollen metabarcoding is a valuable tool
 to reveal the impact of pollinator foraging behavior on plant-pollinator
 interaction networks in fragmented landscapes.
description: <strong>Data collection</strong>
 Pollinators were collected during the summer of 2021 in two
 sampling periods: 24 June – 9 July and 2 – 13 August. These periods were
 chosen to coincide with the two main flowering peaks in the study area.
 Each of the nineteen dune slacks (Fig. S1) was sampled once during each
 period. Sampling took place between 10:00 and 17:00 under sunny weather
 conditions, with air temperature above 16°C, wind speeds below 3 Bft, and
 cloud cover below 90%, using manual netting, with two persons actively
 walking through the entire dune slacks for two hours. Every insect species
 that was observed foraging on flowers was considered a pollinator. In
 total 3307 pollinators were collected. After capturing, pollinators were
 immediately placed in individual tubes, filled with 100% ethanol to kill
 and preserve them. Whenever possible, pollinators were identified to
 species level using sequential identification keys. If species
 identification was not possible, individuals were assigned to the highest
 possible taxon (85% of samples identified to species level). The remaining
 15% of individuals were mostly identified to genus level, except that
 dipterans of the section Schizophora were identified only to family level
 and further classified into morphotypes. Although all flying insects were
 collected in the field, we often captured only one or two individuals per
 site for rare species. To avoid over-representing common species in the
 pollen analysis and to ensure comparability across taxa, we standardized
 our sampling effort by limiting the analysis to three individuals for
 common species (see section ‘DNA extraction and sequencing’
 below). During both sampling periods, the plant
 community was characterized using 1 m² quadrants, in which the number of
 open flowers for each flowering plant species was recorded. We set a
 minimum of 3 quadrats per dune slack, and additional quadrats were added
 until no new plant species with a relative floral abundance greater than
 10% was encountered. As a result, the number of quadrats varied between
 three and twelve, depending on the size of the dune slack. Floral
 diversity was calculated as the exponent of the Shannon entropy:
                                                       
                              
 <sup>1</sup>D=exp⁡(-∑<sub>i=1</sub><sup>S</sup> p<sub>i</sub>  ln⁡ p<sub>i</sub>),                                                     (1) where <em>S</em> is the total number of flowering species and <em>p<sub>i</sub></em> the proportional abundance of species <em>i</em>. For species in the families Asteraceae and Apiaceae, the number of pseudanthia or umbels was counted instead of the number of flowers to assess their relative abundance. <strong>Patch characteristics</strong> For each dune slack, patch fragmentation was quantified using two metrics: dune slack area (m²) and Hanski’s Connectivity Index (Hanski, 1994), calculated for patch <em>i</em> as , where <em>d<sub>ij</sub></em> is the edge-to-edge distance between all pairs of patches <em>i</em> and <em>j</em>, and <em>A<sub>j</sub></em> is the patch area of patch <em>j</em> in m². For the constant <em>a</em>,** a value of 2.5, which assumes an average foraging distance of 400 m, which is a realistic distance for many pollinator species (Zurbuchen et al., 2010, Desaegher et al. 2022), and was chosen based on Devriese et al. (2024a), who tested other distances as well. The parameter <em>b</em> was set at 0.5, based on Moilanen & Nieminen (2002). Both measures were calculated using the R package “metacapa” (Strimas-Mackey & Brodie, 2018).   <strong>DNA extraction and sequencing</strong> Pollen was extracted from the bodies of the pollinators by manually shaking the tubes in which they were stored. After removing the insects for further identification, the ethanol containing pollen was pooled according to pollinator species for each site and sampling period. In order to prevent uneven sampling distorting the perceived dietary breadth and the number of hetero-specific pollen transported, no more than three individuals per insect species were randomly selected for each site and sampling period for common species. For rarer species, one or two individuals were included as collected. This resulted in 1055 samples. Afterwards, the ethanol was evaporated using a concentrator, leaving behind a dry pollen pellet. DNA was extracted from the pollen using the Plant/Fungi DNA Isolation 96 Well kit (Norgen Biotek). First, the pellets were resuspended in 200 µL of Lysis Buffer L and transferred to a 96-well plate with ceramic beads. Samples were homogenized using a bead mill. After mechanical homogenization, an additional 200 µL of Lysis Buffer L was added to each well, and DNA extraction was performed following the protocol supplied by the manufacturer. Library preparation and sequencing were performed as described by Devriese et al. (2024b), using the primer pair ITS2-S2F/ITS4R (White et al., 1990; Chen et al., 2010). Sequencing was performed at Genomics Core Leuven using the Illumina MiSeq platform with the v2 500-cycle reagent kit (Illumina, San Diego, CA, USA). <strong>Bioinformatics</strong> The filtering and construction of the zOTU table were performed using the pipeline developed by Leonhardt et al. (2022). Forward and reverse paired-end reads were merged with a minimal overlap of 20 base pairs and a maximal difference of 10 base pairs. Only merged sequences between 200 and 400 base pairs were retained, with a maximum expected error rate of one. Denoising was performed with the α parameter set to two (Edgar, 2016). Zero‐radius operational taxonomic units (zOTUs) were identified using the unoise3 algorithm with the default settings from USEARCH (Edgar, 2013). Taxonomy was assigned at a 97% similarity threshold, using the ITS2 European Countries database (Quaresma et al., 2024) for Belgium, the Netherlands, France, Germany, and a global database. Unassigned zOTUs were manually assigned using a BLAST search in NCBI gene bank (https://blast.ncbi.nlm.nih.gov/Blast.cgi). To remove false reads, zOTUs with fewer than 16 reads or a relative read abundance below 0.1%  in a sample were removed using the phyloseq package in R version 4.0.3 (McMurdie & Holmes, 2013). Finally, to filter out tag jumping, sequences contributing less than 0.01% of the total reads of a plant species were removed from the dataset (Drake et al., 2022; Rodriguez-Martinez et al., 2023). All observed species were validated against the database “Waarnemingen.be” to confirm their presence in the study region. Nine species that were not found in the database, were removed. <strong>Network properties</strong> For each of the 19 sampled dune slacks, based on the metabarcoding results, a binary interaction network was constructed at species-species level, by identifying all links between plants and insects, combing both sampling periods. Connectance (<em>C</em>), i.e., the proportion of all possible links that are realized, and nestedness based on overlap and decreasing fill (NODF), which allows to compare networks of different sizes (Almeida-Neto et al., 2008) were calculated using the “networklevel” function from the bipartite package (Dormann et al., 2022). Modularity <em>Q</em> (Newman, 2006), was calculated using the simulated annealing algorithm from rnetcarto (Doulcier & Stouffer, 2023), based on Guimera et al. (2007). To assess hierarchical compound network topology, within-module nestedness (NODFsm) and between-module nestedness (NODFdm) were calculated following Flores et al. (2013) and Felix et al. (2017). Here, instead of averaging the nestedness score between pairs of species over the whole network, as is the case for NODF, the nestedness scores are averaged separately, for pairs of species from the same module (NODFsm), and pairs of species from different modules (NODFdm). All metrics were calculated for binary networks and tested for spatial autocorrelation with Mantel correlograms using the function “mantel.correlog “ from the R package “vegan” (Table S1) (Oksanen et al. 2022). To identify patterns in module compositions, cluster analysis was performed. A matrix was constructed indicating the module assignment of each species across the networks. For plant species, only those that occurred in at least 16 out of the 19 networks were included, while for pollinators, all species that occurred in at least 8 networks were included. From these matrices, distance matrices were created using the Gower distance in the R package “cluster” (R Core Team, 2022). Subsequently, clustering was performed using Ward’s hierarchical agglomerative clustering method (Murtagh & Legendre, 2014), as implemented in the function hclust of the ”stats” package in R (R Core Team, 2022). Furthermore, to identify species and species groups that connect modules to each other, among-module connectivity <em>c</em> was calculated for each species in each network, using the function “netcarto” from the package “rnetcarto” (Doulcier & Stouffer, 2023). Among module connectivity for species <em>i</em>, also known as the participation coefficient, is calculated as: c<sub>i=1</sub>- ∑<sub>s=1</sub><sup>N</sup>(k<sub>is</sub>/k<sub>i</sub>)<sup>2</sup>  ,  where <em>N</em> is the number of modules, <em>k<sub>is</sub></em> the number of interactions within module <em>s</em> and <em>k<sub>i</sub></em>  the total number of interactions for species <em>i</em> (Guimerà & Amaral, 2005). To identify species that define module composition, the within-module degree <em>z</em> was calculated as: z<sub>i</sub>=(k<sub>i</sub>-k ̅<sub>si</sub>)/σk<sub>si</sub>  , with <em>k<sub>i</sub></em> the number of interactions of species <em>i</em> within module <em>si</em>,  the mean number of interactions for all species in module <em>si</em> and**  the standard deviation of <em>k</em> in module <em>si</em> (Guimerà & Amaral, 2005). <strong>Null models</strong> To account for the effects of network size and connectance on other network metrics and to standardize these metrics across networks, z-scores were calculated (Ulrich et al., 2009). Different null models were used to calculate z-scores as:   ,                                             Z=(x<sub>o</sub>-x<sub>null</sub> ) ̅/σ<sub>null</sub>                              (2) where <em>x</em><sub>0</sub> represents the observed network parameter, <em>x</em><sub>null</sub> is the mean network parameter from the null models, and s<sub>null</sub> is the standard deviation from the null models. In addition, this approach allowed to test whether the observed networks were ecologically relevant or simply the result of random processes. As not all network metrics have the same conditions, distinct null models, with different constraints, were used. For modularity and nestedness, the “Curveball” algorithm from Strona et al. (2014) was used, with 10000 iterations, 1000 burn-in steps, and a thinning factor of 500. This null model is a binary sequential swap algorithm that preserves connectance, as well as the degree sequences of both plants and pollinators. For compound structure, NODFsm and NODFdm, the Curveball algorithm is not applicable and therefore the proportional-restricted null models of Felix et al. (2017) were used with 10000 iterations. These null models preserve the connectance and modular structure of the original network while keeping the degree distribution proportional to that of the original network (Felix et al. 2022). The obtained *p-*values were corrected using the Benjamini-Hochberg adjustment (Benjamini & Hochberg, 1995). <strong>Effects of fragmentation</strong> These analyses resulted in six different network metrics: connectance (C), nestedness (NODF), modularity <em>Q</em> , within-module nestedness (NODFsm), between-module nestedness (NODFdm) and the ratio between within-module and between-module nestedness (NODFsm/NODFdm). In addition, for each measure, relative measures were calculated using the above mentioned null models: relative modularity (zM), relative nestedness (zNODF), relative within-module nestedness (zNODFsm), and relative between-module nestedness (zNODFdm). To assess whether topological features of networks were significantly related to patch characteristics, linear models were constructed relating each network metric to log-transformed patch area, Hanski’s connectivity index and floral diversity as independent variables. Linear models were analyzed using the R package “stats” and “car” (R Core Team, 2022). All models met the assumptions of linearity, homoscedasticity and normality of the residuals. For NODFsm and NODFdm, along with their respective <em>z</em>-scores and ratios, only networks that were significantly more modular than predicted by the null models were included in the analysis.
description: # Interactions outside local patches contribute to the compound topology
 of plant-pollinator networks in fragmented dune slacks
 [https://doi.org/10.5061/dryad.sj3tx96dp](https://doi.org/10.5061/dryad.sj3tx96dp) ## Description of the data and file structure **Interactions outside local patches preserve compound topology of plant-pollinator networks in fragmented dune slacks** In this study, we used DNA metabarcoding to construct pollen transport networks across nineteen dune slacks along a fragmentation gradient at the Belgian coast and  investigated the impact of patch area, connectivity, and floral diversity on network size, connectance, specialization, nestedness and modularity. Networks showed a low degree of specialization and none of them were significantly nested, while most were significantly modular with a hierarchical compound structure. Importantly, no significant relationships were found between network measures and patch area or connectivity, indicating that local patch characteristics had little impact on network structure. We conclude that plant-pollinator networks in fragmented landscapes with a permeable matrix are able to maintain network characteristics trough interactions with plants outside the local patch, and that pollen metabarcoding is a valuable tool to reveal the impact of pollinator foraging behaviour on plant-pollinator interaction networks in fragmented landscapes. ### Files and variables #### File: Raw_sequences.zip **Description:** The file consists of three separate folders which each represent a separate sequencing library. They all contain raw sequencing files, fastq files are named according to the sampling code of "Extraction_samples.csv" (see further below). Some subfolders have the files named without the preceding "A". ## Code/software `Filtering.zip` contains the scripts and necessary data for the filtering of the output of the bioinformatics. It is located on Zenodo: * raw ASV table "asv_table.merge.txt", * taxonomy file: "taxonomy.txt",and * sample information: "Extraction_samples.csv" * an R-script to filter the data: datafiltering.R `processing_and_plotting_(3).zip` contains two r scripts for all statistical analysis and plotting of the data "dataprocessingcode.R" and "netmet.R", which is a dependable for the first script. * Extraction_samplesf.csv, is a condensed sample information file, where failed samples have been filtered out, * Flowc.scv contians floral counts per patch, * JA_publ, contains patch information (site, sampling month area (m^2) and hanski connectivity) , * Pollinatormatrix.csv contains the pollinator counts per patch, * Result_matirx_extr.csv is a sample-plant species matrix, retrieved after the filtering, * speclas.csv contains weather plant species are classified as occuring in dune slacks (s), the landscape matrix (m) or gardens (t) First run the "netmet" script before running "dataprocessingcode". To run "netmet you need dependancies that can be retrieved from: [https://github.com/gabrielmfelix/Restricted-Null-Model](https://github.com/gabrielmfelix/Restricted-Null-Model) Some translations for non Duch speaking users: "Pan" = Dune slack, eg patch ID "Maand"= Month, eg sampling period, J=July, A=August "Soort" = Species