A network of 112 bomb crater ponds forming a well-delineated pondscape relatively isolated from other waterbodies is situated in the Kiskunság area, Central Hungary (Fig. 1). These ponds were likely created during World War II by mistargeted bombing on a sodic meadow. They are saline waters mostly dominated by sodium carbonates and hydrocarbonates, and despite being the same age and in close proximity to each other, they vary in many of their environmental and morphological characteristics including hydroperiod, surface area, depth and vegetation cover (Vadet al. , 2017). The ponds are not physically connected thus dispersal is expected to occur via wind, animal vectors or active movement of the organisms.Figure 1. A) Location of the study region near Apaj in Hungary. B) Drone photo of a smaller part of the bomb crater pond network. C) Location of the 54 bomb crater ponds with their closeness centrality scores indicated by a colour gradient (Google Earth Pro, 2014)

Fieldwork was carried out between 7 and 9 May 2014. We surveyed all ponds holding water (i.e., omitting the smallest and already dry ephemeral habitats), resulting in a total of 54 ponds. The average distance between these ponds is 285 m and all habitats are situated within a radius of 0.5 km (min. and max. distances between the centre point of ponds: 9–863 m). A range of physical and chemical parameters was recorded on site including water depth (cm), diameter (m), and open surface area (%), emergent and submerged vegetation (%). Conductivity (μS cm⁻¹) and pH were measured using a field multimeter (Eutech CyberScan PCD 650). Per pond, a total of 10 L of water was collected and mixed from 10 randomly chosen points, 1 L of this composite water sample was taken for water chemistry analysis and community sequencing of pro- and microeukaryotes. 1-50 mL (depending on turbidity) of this water sample was filtered through a nitrocellulose membrane filter (Ø 47 mm, 0.22 μm pore size). These filters were stored at −20 °C until DNA extraction. Concentrations of total phosphorus (TP, mg L^-1), total suspended solids (TSS, mg L^-1), dissolved inorganic nitrogen (DIN, mg L^-1), calcium (Ca²⁺, mg L^-1), and chlorophyll a (Chl-a ćg L^-1) were measured in the lab as detailed in Vadet al. , (2017).

To collect zooplankton samples, 10 L of water was randomly collected from the open water and sieved through a 45-μm mesh. Macroinvertebrates were sampled using sweep netting standardised to 3 minutes (frame-size: 0.25x0.25 m², mesh size: 500μm), for which all microhabitats present in a pond were included. Zooplankton and macroinvertebrate samples were preserved in 70% ethanol for further analysis. For the amphibian survey, hand netting maximised to 15 minutes, visual searches, and dip netting for tadpoles and newt larvae were applied. All identified specimens were released back to their habitat. For a more detailed description of the study area and sampling procedures, see Vad et al. (2017)

The isolation of community DNA from the membrane filters was carried out using the PowerSoil® DNA Isolation Kit (MO BIO Laboratories Inc.). Prokaryotic 16S and microeukaryotic 18S rRNA gene amplification and sequencing were performed by LGC Genomics (Berlin, Germany). For gene amplification, the following primer pairs were applied: EMBf 515F (GTGYCAGCMGCCGCGGTAA, Parada et al., 2016) – EMBr 806R (GGACTACNVGGGTWTCTAAT, Apprill et al., 2015) for the V4 region of the prokaryotic 16S rRNA gene and UnivF-1183mod (AATTTGACTCAACRCGGG) – UnivR-1443mod (GRGCATCACAGACCTG) (Ray et al., 2016) for the V7 region of the eukaryotic 18S rRNA gene. Sequencing was carried out on an Illumina MiSeq platform. For a detailed description of the entire procedure, see Szabó et al. (2022). Amplicon readsets were analysed using mothur v1.43.0 (Schloss et al. , 2009) involving, sequence processing, taxonomic assignments, and OTU picking with the MiSeq SOP as a reference (http://www.mothur.org/wiki/MiSeq_SOP; Kozich et al. , 2013, downloaded on 12^thNovember 2020). Additional quality filtering steps were implemented to minimize the presence of sequence artefacts. These steps included the adjustment of the deltaq parameter to 10 in the ’make.contigs’ command, primer removal from both ends of the sequences, chimera identification and removal using the mothur-implemented version of VSEARCH and the exclusion of singleton reads. De-noising was carried out using mothur’s ‘pre.cluster’ command with the suggested 2 bp difference cutoff. Read alignment and taxonomic assignment were performed using the ARB-SILVA SSU Ref NR 138 reference database with a minimum bootstrap confidence score of 80 (Quast et al. , 2012). Non-primer-specific taxonomic groups (‘Chloroplast’, ‘Mitochondria’ and ‘unknown’) were excluded from the dataset.

OTUs were selected at the 99% similarity threshold. Taxonomic assignment of 18S rRNA gene OTUs was performed using the ‘classify.otu’ command with the PR2 v4.12.0 reference database (Guillou et al. , 2012). 18S rRNA gene OTUs assigned to taxa Streptophyta, Metazoa, Ascomycota, and Basidiomycota were excluded from the microeukaryotic dataset. For the taxonomic assignment of prokaryotic OTUs, the TaxAss software (Rohwer et al. , 2018) was used with default parameters and the FreshTrain (15 June 2020 release) and ARB-SILVA SSU Ref NR 138 databases. Subsequently, both 16S and 18S OTU sets were rarefied to the read number of the sample having the lowest sequence count. Samples BC40 and BC105 were excluded from the eukaryotic dataset due to a low read count after filtering.

Zooplankton was identified to the lowest possible taxonomic level by microscopic analysis (generally to species), except for bdelloid rotifers which were treated as a single group. Macroinvertebrates were also identified to the lowest possible taxonomic level (generally to family, genus and species levels, depending on the higher taxa) using relevant identification keys. As passively dispersing macroinvertebrates were represented only by four taxa (Vad et al. , 2017), they were excluded from the subsequent analyses. For a detailed description of the sample processing and the list of identified taxa, see Vad et al. , (2017).

We focused on the following organism groups in our study (Table S1; Supporting information): prokaryotes, phototrophic microeukaryotes, heterotrophic microeukaryotes, crustacean zooplankton (copepods and cladocerans), rotifer zooplankton, dipterans, other macroinvertebrates, and amphibians-reptilians. For community matrices of prokaryotes and microeukaryotes, we used the rarified 16S and 18S OTU tables. The 18S dataset was split into phototrophs and heterotrophs based on phyla using the ‘phyloseq’ package v1.36.0 (McMurdie & Holmes, 2013) in R 4.1.0 (R Core Team, 2021). Phototrophs mostly included groups capable of photosynthesis, while heterotrophs consisted of the remaining groups. Dipterans were treated separately from other macroinvertebrates as they are considered weak active dispersers (Bilton, Freeland, & Okamura, 2001; Heino, 2013). Other macroinvertebrates included all actively dispersing groups considered intermediate or strong aerial dispersers with flying adults according to Heino (2013).

In order to assess the relative importance of environmental and spatial predictors in explaining patterns of taxonomic richness and community similarity, we carried out three sets of variance partitioning analyses using the varpart function of the ‘vegan’ package (Oksanenet al. , 2020), where either network position (in the case of taxonomic richness) or spatial distance (in the case of taxonomic richness and community similarity) was included as the spatial predictor. We used data from all 54 ponds in all three cases, except for phototrophic and heterotrophic microeukaryotes for which data was available from 52 ponds only. In the environmental dataset, we excluded highly correlated (i.e., Pearson’s r>0.75) environmental variables (turbidity, percentage of reed cover) prior to the analyses. Secchi depth and submerged vegetation cover were also excluded, as Secchi depth delivers the same information as TSS but on a less precise scale, while the percentage of submerged vegetation cover almost exclusively contained zeros. Our final environmental dataset thus contained 10 variables, which were tested for normality and those deviating from it (Shapiro-Wilk test p <0.05) were transformed (conductivity, pH, depth - untransformed; DIN, TP, TSS - natural log; Ca²⁺, surface area - square root; chlorophyll-a- log(x+1); open surface area - cube transformed).

First, we analysed how taxonomic richness is predicted by the local environment and the relative position of ponds within the habitat network. For this, partial multiple linear regressions for each studied organism group were carried out using the rda function of the ‘vegan’ package (Oksanen et al. , 2020). As a response variable, richness was included either as OTU (pro- and microeukaryotes) or taxonomic richness (other groups), but hereinafter, we refer to these as ‘taxonomic richness’ for simplicity. Taxonomic richness data of crustacean zooplankton and amphibians-reptilians were square-root transformed to normalise model residuals (inspected via diagnostic plots), while untransformed values were used in the other organism groups. As environmental predictors, we used the first two axes of a Principal Component Analysis (PCA) constructed for the environmental variables, following z-score standardisation. For all subsequent analyses, site scores of PC1 and PC2 axes (explaining >65% of the total variation, Fig. S1, Supporting information) were extracted and used as environmental data. As a spatial predictor, we calculated the closeness centrality index of each pond (Fig. 1 C) based on the spatial distance matrix of all ponds using the ‘igraph’, ‘fields’ and ‘reshape2’ packages (Csárdi & Nepusz, 2006; Wickham, 2007; Nychkaet al. , 2021). To test for the possible significance of the unique effects of environment and centrality, a permutation test (2000 permutations) was run for each partial multiple linear regression model using the function anova .

We run a second set of variance partitioning (partial multiple linear regression models) for each organism group separately to test if and how taxonomic richness is affected by the environmental variables (PC1 and PC2 axes) and the number of neighbouring ponds. In these models, environmental variables were included as one set of predictors the same way as described above, and the number of ponds within a radius with increasing distance from each local pond over a scale of 0-800 m, with 10 m increments was involved as the spatial predictor. The unique effects of space and environment were calculated the same way as in the first set of variance partitioning 80 times in total for each organism group. The unique effect of space along threshold distance (cut-off distance for calculating pond numbers) was plotted with the help of General Additive Models (GAMs, formula = y ~ s(x, k =5)) for each taxonomic group separately with the package ‘ggplot2’ (Wickham, 2016).

Third, we tested how metacommunity structure is predicted by the local environment and spatial configuration. Since we did not have abundance data for all organism groups, statistical analyses were run on presence-absence community data for comparability. Taxa occurring in less than three ponds were omitted given their minor contribution to the overall similarity in the community dataset. Then, any pond which contained no taxa was removed. Dissimilarity matrices were calculated based on Sørensen dissimilarity with the vegdist function of the ‘vegan’ package (Oksanen et al. , 2020). Here, the initial pool of environmental predictors was represented by the 10 transformed environmental predictors, while the spatial predictors were Moran’s Eigenvector Maps (MEM). These were constructed based on the spatial distance matrix using the dbMEMfunction from the ‘adespatial’ package (Dray et al. , 2022) and only the 15 positive MEMs were retained (Bauman et al. , 2018; Borcard, Gillet, & Legendre, 2018; Fig. S2; Supporting information). To select significant variables of each set of explanatory variables (environment and space), we applied stepwise selection both for the transformed environmental variables and for the positive MEM eigenvectors using the ordistep function (direction=both, perm.max=2000) of the ‘vegan’ package (Oksanen et al. , 2020). These selection steps were done separately for each organism group resulting in different sets of retained environmental variables and MEMs (Table S2, Supporting information). We carried out a set of partial distance-based redundancy analyses (dbRDA) using the capscalefunction of the ‘vegan’ package (Oksanen et al. , 2020) followed by significance testing with a permutation test (2000 permutations) using the function anova .

All statistical analyses were carried out in R 4.1.0 (R Core Team, 2021) and figures were created using the ‘ggplot2’ (Wickham, 2016), ‘pubbr’ (Kassambara, 2020) and ‘car’ packages (Fox & Weisberg, 2019).