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.

Data collection

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:                             

                                                        ¹D=exp⁡(-∑_i=1^S p_i  ln⁡ p_i),                                                     (1)

where S is the total number of flowering species and p_i the proportional abundance of species i. 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.

Patch characteristics

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 i as , where d_ij is the edge-to-edge distance between all pairs of patches i and j, and A_j is the patch area of patch j in m². For the constant a,** 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 b was set at 0.5, based on Moilanen & Nieminen (2002). Both measures were calculated using the R package “metacapa” (Strimas-Mackey & Brodie, 2018).

DNA extraction and sequencing

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).

Bioinformatics

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.

Network properties

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 (C), 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 Q (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 c was calculated for each species in each network, using the function “netcarto” from the package “rnetcarto” (Doulcier & Stouffer, 2023). Among module connectivity for species i, also known as the participation coefficient, is calculated as: c_i=1- ∑_s=1^N(k_is/k_i)²  ,  where N is the number of modules, k_is the number of interactions within module s and k_i  the total number of interactions for species i (Guimerà & Amaral, 2005). To identify species that define module composition, the within-module degree z was calculated as: z_i=(k_i-k ̅_si)/σk_si  , with k_i the number of interactions of species i within module si,  the mean number of interactions for all species in module si and**  the standard deviation of k in module si (Guimerà & Amaral, 2005).

Null models

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_o-x_null ) ̅/σ_null                              (2)

where x₀ represents the observed network parameter, x_null is the mean network parameter from the null models, and s_null 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).

Effects of fragmentation

These analyses resulted in six different network metrics: connectance (C), nestedness (NODF), modularity Q , 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 z-scores and ratios, only networks that were significantly more modular than predicted by the null models were included in the analysis.