Linking biodiversity and ecological function through extensive microeukaryotic movement across different habitats in six urban parks

Highly diverse but divergent microeukaryotes dwell in all types of habitats in urban park ecosystems. Extensive microbial migration occurs between both terrestrial and aquatic habitats. Microbial movement is beneficial to the maintenance of biodiversity and the exchange of functional guilds.


Sample processing
Six urban parks, including Shuanglongtan park (0.55 square kilometer), Dapingshan park (1.50 square kilometer), Huli park (0.11 square kilometer), Shangli park (0.48 square kilometer), Dalunshan park (0.36 square kilometer), and Xiangshan park (10.38 square kilometer) were sampled in July 2020.The moss, surface soil (top 5 cm), and surface sediment (5 cm) were scooped up by a trowel.The litter or detritus in the tree hole, including dust and the small tree bark, were collected.The sample represented 100 g of substrates or 2.5 L of surface waters (0.5 m), which were preserved in ice boxes before treatment in the laboratory.The water samples (500 mL) were pre-filtered through 200 µm mesh to remove larger particles and metazoans and subsequently filtered through a 0.22 µm pore-size polycarbonate filter (47 mm diameter, Millipore, Billerica, MA, USA). 2 g of moss was suspended in 40 mL phosphate buffer saline (PBS), and shaken for 2 h under 180 rpm at 30 ℃.The solution was treated with ultrasound for 10 minutes, then left to stand for 2 h.After that, microbes were filtered through a 0.22 mm pore-size polycarbonate filter.All samples were stored at -80 °C until DNA extraction.

DNA extraction and sequencing
For sediments, soils, and tree holes, 0.5 g samples were used to extract DNA.For moss and water samples, filters were cut into small pieces with a sterilized cutter.The total DNA was extracted using the Fast DNA SPIN Kit and the Fast Prep Instrument (MP Biomedicals, Solon, OH, USA) following the manufacturer's instructions.For the whole microeukaryotic community, the V4 or V9 region of the eukaryotic 18S rRNA gene was amplified using the primer pair 547F/967R [1] or 1380F/1510R [2], respectively.A 30 μL polymerase chain reaction (PCR) reaction included 15 μL of Phusion High-Fidelity PCR Master Mix (New England Biolabs, Beverly, MA, USA), 0.2 μM of each primer, and 10 ng of DNA.For the V4 region, the PCR reaction comprised of an initial denaturing step at 95°C for 5 min, followed by 30 cycles of 94°C for 30 s, 45°C for 45 s, and 72°C for 60 s.For the V9 region, the PCR protocol employed an initial denaturing step at 95°C for 5 min, followed by 30 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s.At the end of the amplification, the amplicons were subjected to a final 10-min extension at 72 °C.The PCR reaction was performed in triplicate per sample, and then PCR products from triplicate reactions were pooled and gel-purified.All libraries were sequenced on the Illumina platform (Illumina Inc., San Diego, CA, USA) using a paired-end (2 × 250 bp for V4 region; 2 × 150 bp for V9 region) strategy.

Sequence quality control
Paired-end reads were assigned to each sample based on their unique barcodes.Quality filtering was performed by eliminating the low-quality reads with QIIME V1.9.1 [3].
Default settings for quality control processing were used: the maximum number of consecutive low-quality base calls allowed before truncating a read = 3, the minimum number of consecutive high-quality base calls to retain read = 0.75; last quality score considered low quality = 3; the maximum number of ambiguous characters allowed in a sequence = 0.The chimera sequences were detected and removed by using the UCHIME de novo algorithm [4].

Statistical analyses
All samples from V4 and V9 regions were merged at the supergroup level, and hierarchical clustering using the "ward.D2" method based on Bray-Curtis dissimilarity was implemented using R software (version 4.1.0)[5] and visualized by the "ggtree" package [6].Permutational multivariate analysis of variance (PERMANOVA) at the supergroup level was conducted to explore the relative contribution of amplicon region, park, and habitat to the difference at the community level.Rarefaction curves were calculated based on individual samples and each of the habitats to explore whether the sequencing depth was sufficient to cover the majority of microeukaryotic taxa.Indices of α-diversity, including richness and phylogenetic diversity (PD) were calculated.
Niche breadths of microeukaryotic groups were estimated through the Levins' niche breadth index by the "spaa" package [7].Ecological niche breadth of a species represents its ability to utilize a range of resources, and a wider niche breadth is generally considered to have higher metabolic flexibility [8].Non-metric multidimensional scaling (NMDS) ordination and analysis of similarities (ANOSIM) were used to investigate differences in community compositions among different groups (five habitats and six parks).
The community compositions between samples were analyzed based on the Bray-Curtis dissimilarity.These analyses were performed by the "vegan" [9], "picante" [10], and "ape" [11] packages.Spearman correlation of Bray-Curtis dissimilarity was used to test the significance of the correlation between V4 and V9-based community compositions at different habitats.The overall β-diversity between any two habitats was partitioned into turnover and nestedness components by the "betapart" package [12].
The significance test between different habitats or sequencing regions (V4 vs V9) was performed by the Wilcoxon test.
We evaluated the fit of the Sloan neutral community model for microeukaryotic community compositions to determine the potential importance of neutral processes on community assembly [13].The value Nm determines the correlation between occurrence frequency and regional relative abundance, with N indicates the community size and m being the immigration rate.The R 2 represents the overall fit to the neutral model.Models were run following a published R script [14,15].To further infer the ecological assembly processes, we calculated the β-nearest taxon index (βNTI) [16].βNTI > 2 indicates variable selection of deterministic process; βNTI < −2 indicates homogeneous selection of the deterministic process.If the |βNTI | < 2, the stochastic processes are deemed to drive observed differences.Then, community assembly was further distinguished using the modified Raup-Crick index (RCBray), which is able to disentangle dispersal limitation (indicated by RCBray > 0.95) and homogenizing dispersal (indicated by RCBray < −0.95).If the |RCBray| < 0.95, an "undominated" fraction, which contains weak selection, weak dispersal, diversification, and drift shaped community structure.FEAST analysis uses an expectation-maximization algorithm and can track the contribution of numerous potential source environments simultaneously [17].Here, the five habitats in pairs from six urban parks were selected as unknown sinks and potential sources, respectively.There were 18 replicates for each habitat (3 replicates for each habitat in one park × 6 parks), hence, there were 18 replicates for each sink-sources matrix.The input file was prepared under the guidance of the FEAST package, and 1000 iterations were used to perform FEAST analysis.A, B,     E) and V9 (C, D, F) regions, respectively.Wilcoxon test was performed to verify whether a zOTU's abundance significantly differs between each pair of habitats.We assume that species may have migrated between the two habitats if there is no significant difference in the abundance of the same species in these two habitats.

FigureFigure S2 Figure S3 Figure S4
Figure S1 Rarefaction curves of microeukaryotic community based on V4 and V9 regions of 18S rRNA gene, respectively.(A, B) The individual samples; (C, D) The combined sets of five-habitat samples.

Figure S5 FitFigure S6 Figure S7 Figure S8
Figure S5 Fit of the neutral community model of microeukaryotic community from urban parks.The predicted occurrence frequencies along the different habitats.Single (A) and multi-habitat (B) based on V4 region-communities; Single (C) and multi-habitat (D) based on V9 region-communities.R 2 indicates the fit to the neutral model and is shown in the bar plot.Negative R 2 values can occur when there is no fit to the model.Nm indicates the community size time immigration and is shown in a line chart.
(A, C) and relative percentages (B, D) of migrated zOTUs between pairs of habitats.(E, F) Total relative abundance of migrated zOTUs in each pair of habitats.The x-axis is each habitat pair.

Figure S10
Figure S10 Potential migrated fungal taxa with assigned trophic mode from 18S rRNA gene V4 (A, B) and V9 (C, D) regions, respectively.Numbers (A, C) and relative percentages (B, D) of migrated zOTUs between pairs of habitats.

Table S2
Observed relative abundances of testate amoeba cells from microscopic analysis, and sequence numbers and relative abundances of testate amoebae from amplifying 18S rRNA gene V4 and V9 regions.Pearson and Spearman correlations of relative abundance are calculated between microscopic analysis and V4-based testate