Heterogeneity across Neotropical aquatic environments affects prokaryotic and eukaryotic biodiversity based on environmental DNA

Characterizing biological communities and knowledge on the distribution of biodiversity allows the assessment of ecological quality. This provides valuable information for conservation biology and monitoring purposes. While obtaining such data has been challenging in the past, environmental DNA (eDNA) sampling represents a promising tool to describe biodiversity on a broad taxonomic scale. In this study, we provide the first broad-scale biodiversity assessment for ten Neotropical water bodies in Nicaragua (a major river, two great lakes

On a global scale, the distribution of biodiversity can be affected by multiple abiotic factors such as latitude and altitude, and also precipitation, for example (reviewed in Gaston, 2000).Since such patterns are often consistent for different taxonomic groups (Peters et al., 2016;Stevens, 1989), one could expect positive covariance in species richness across different taxa.While positive associations have been shown, for example, for different groups of Australian vertebrates (Pianka & Schall, 1981) or across trophic levels in aquatic habitats (Andersen et al., 2020), concordance in species richness can be highly variable depending on the investigated taxa, their trophic level (Mjelde et al., 2012), and geographic scale (Flather et al., 1997;van Jaarsveld et al., 1998;Prendergast et al., 1993;Westgate et al., 2014).However, due to the limited capacity of effectively quantifying species richness and diversity across different taxonomic groups and ecosystems, it has been challenging to reach a general consensus on the factors that affect the distribution of biodiversity.
Freshwater ecosystems generally show high levels of species richness compared with oceanic and terrestrial ecosystems; while they cover only about 2.3% of our planet's surface, they harbor around 5%-6% of all described species (and almost 9.5% of all described animal species; Dudgeon et al., 2006;Grosberg et al., 2012;Reid et al., 2019).Species richness of freshwater fishes and amphibians tends to increase toward lower latitudes, and particularly, Central and South America are hot spots for freshwater biodiversity (Abell et al., 2008;Pyron & Wiens, 2013).Freshwater ecosystems are also highly threatened by the loss of biodiversity due to a wide range of stressors (Reid et al., 2019;Su et al., 2021), especially in tropical regions (Huete-Pérez et al., 2016).For example, the construction of an interoceanic shipping canal was feared to imperil the freshwater biodiversity of several water bodies in Nicaragua (Härer et al., 2017;Huete-Pérez et al., 2015;Huete-Pérez et al., 2013, 2016).This emphasizes the need to better understand the distribution of species richness and diversity across different environments in order to develop efficient conservation strategies.
Surveying biodiversity at a broad taxonomic scale by traditional methods poses numerous technical challenges, for example, invasive sampling and microscopic identification (Ficetola et al., 2008).
These challenges can be overcome by environmental DNA (eDNA) sampling (Deiner et al., 2017;Huerlimann et al., 2020;Nguyen et al., 2020;Taberlet et al., 2012).But, eDNA metabarcoding comes with challenges if there are limited data for reference sequences, which might hinder identification on the species level (Jackman et al., 2021).Nevertheless, metabarcoding of eDNA provides a non-invasive and highly efficient approach to detect a wide range of prokaryotic and eukaryotic taxa and it is being established as an important tool for biodiversity monitoring, particularly in aquatic environments (Deiner et al., 2016(Deiner et al., , 2017;;Pawlowski et al., 2020;Rees et al., 2014).However, more biodiversity monitoring using eDNA sampling is conducted in temperate environments than in tropical environments (Huerlimann et al., 2020).This is because the integrity of eDNA is thought to be unstable due to the higher impact of environmental stressors (e.g., rainfall, temperature, and UV level).It has just recently been suggested that eDNA from tropical areas is stable enough in the environment for biodiversity monitoring, but only few empirical examples have been published to date (Dommain et al., 2020;Robson et al., 2016).Therefore, more studies are needed to further explore patterns of biodiversity in tropical environments.

| Freshwater environments of Nicaragua
Nicaragua is characterized by a diverse aquatic landscape and is referred to as "the land of lakes and volcanoes."This applies particularly to the western part of the country, which belongs to the Central American Volcanic Arc (Kutterolf et al., 2007).Of particular socio-economic but also scientific interest are the water bodies of the San Juan drainage basin.The two great lakes, Lake Nicaragua and Lake Managua, are rather old with an age of approximately 500,000 to one million years (Bussing, 1976).Lake Nicaragua is the largest freshwater lake in Central America with a surface area of 8264 km 2 ; Lake Managua is much smaller with a surface area of 1042 km².Multiple rivers drain into the two great lakes, and Río San Juan is the only river draining out of Lake Nicaragua at its southeastern end.Río San Juan has a total length of 192 km and connects Lake Nicaragua with the Caribbean Sea.There are also a large number of much smaller (<21.1 km 2 ) and younger (<24,000 years) crater lakes that formed in the calderas of extinguished volcanoes (Figure 1; Barluenga & Meyer, 2010;Kutterolf et al., 2007).Notably, these crater lakes are also much deeper (mean depth: 21.6-142 m) than the great lakes Nicaragua (mean depth: 9 m) and Managua (mean depth: 9.5 m, Table 1).Most crater lakes are surrounded by high rims that completely isolate them from each other and from the great lakes.
Besides their age and size, the lakes also have different chemical and physical characteristics, including turbidity (Torres-Dowdall et al., 2017), pH (7.33-8.69),levels of dissolved oxygen (0.1%-97.3%), and salinity (0.12-4.77ppt; Table 1).These highly variable aquatic environments located within a small geographic range represent a rather unique setting and allow testing whether biodiversity is affected by a range of abiotic factors.Further, this ecological setting is threatened by anthropogenic impacts, such as waste disposal into the lakes or the planned construction of an interoceanic canal, which will also change the physicochemical characteristics of these aquatic environments (Härer et al., 2017;Huete-Pérez et al., 2016;Lacayo et al., 1991).
The influence of abiotic factors on biodiversity remains poorly understood for many organisms across the tree of life (Vilmi et al., 2016).In Nicaraguan lakes and rivers, species diversity has been studied in detail merely for freshwater fish and there is sound knowledge on species assemblages for the lakes investigated in this study (Bussing, 1976;Härer et al., 2017;Torres-Dowdall & Meyer, 2021).Two clear patterns emerge from these data; crater lakes vary considerably in the assemblage of fish they harbor.The two great lakes have a much higher diversity of fish species (more than 65 described species) than crater lakes (2-19 described species; Bussing, 1976;Härer et al., 2017;Torres-Dowdall & Meyer, 2021).The higher diversity in the great lakes is presumably due to the early colonization of these environments by riverine fishes soon after the lakes formed one million years ago (Bussing, 1976).In contrast, colonization of crater lakes by fish is assumed to have occurred more recently from the two great lakes and has been studied comprehensively only for the adaptive radiation of Midas cichlid fish (Amphilophus cf.citrinellus), a model organism for recent speciation in sympatry (Barluenga et al., 2006;Kautt et al., 2016Kautt et al., , 2020) ) and for few other cichlid species (Elmer et al., 2013;Franchini et al., 2017;Xiong et al., 2020).
Colonization events are assumed to have been rare and, to some extent, stochastic.This might be particularly true for large-sized species such as fish.However, little is known about the colonization history and geographic distribution of other organisms, such as invertebrates or microbes.
Here, we utilized eDNA collected from water samples to assess the distribution of biodiversity (measured as α-and β-diversity) across multiple aquatic environments that differ substantially in a range of abiotic factors.These water bodies are located within the San Juan drainage basin of Nicaragua and include Río San Juan, the two great lakes, and seven crater lakes.We specifically tested whether (i) levels of prokaryotic and eukaryotic biodiversity are higher in great lakes compared with crater lakes, as previously described for teleost fishes (Torres-Dowdall & Meyer, 2021), (ii) overall prokaryotic and eukaryotic biodiversity is affected by physicochemical factors that vary across different aquatic environments and (iii) levels of biodiversity are correlated between prokaryotes and eukaryotes but also among phyla within these domains.

| Sampling
Sampling sites included the two great lakes (Nicaragua and Managua), seven volcanic crater lakes (Apoyeque, Apoyo, As.Leon, As.Managua, Masaya, Tiscapa, and Xiloá), and one river (Río San Juan) located in Nicaragua (Figure 1, Table S1).Each water body exhibits unique features in terms of size, depth, age as well as chemical and physical properties (Table 1).It should be noted that we did not collect such information for Río San Juan, which is why it is excluded from some analyses.Four replicates of water samples were collected from the shoreline zone of each water body in 2018, leading to a total of 40 samples.To increase the chances of capturing a more comprehensive estimate of the water body's biological diversity and to avoid bias due to random factors, replicates were taken approximately twenty meters apart.Samples were taken around the same time of day (late mornings) for all water bodies, at intervals of approximately ten minutes.We collected water samples using sterile F I G U R E 1 Geography of the ten water bodies included in this study.Río San Juan flows eastwards out the Atlantic slope of Lake Nicaragua into the Caribbean Sea.All crater lakes are completely isolated water bodies with no inlet or outlet.Crater lakes are indicated by circles, the two great lakes by triangles and Río San Juan by a square symbol.These symbols and colors will be used in subsequent figures

| DNA extraction and amplification
DNA extraction and amplification were carried out under physically separated and specifically dedicated laminar flow hoods to minimize contamination risk.Extractions were done using the QIAGEN DNeasy Blood & Tissue kit (Qiagen) following the manufacturer's instructions.DNA concentrations were measured on a Qubit v2.0 Fluorometer (Thermo Fisher Scientific) and are reported in Table S2.
We amplified the V3-V4 region of the prokaryotic 16S rRNA gene (450 bp; Caporaso et al., 2011) and the V7 region of the eukaryotic 18S rRNA gene (260 bp; Gast et al., 2010;Van de Peer et al., 2000) using well-established primer pairs.DNA amplification was carried out in single reactions, as recommended by Marotz et al. (2019), in two sequential PCRs using the Q5 High-Fidelity polymerase 2× Master Mix (New England Biolabs).After each PCR, the amplified products were purified with HighPrep™ PCR beads (MagBio Genomics).The template DNA was standardized to 2 ng for the first PCR (2 min at 98°C, 10 amplification cycles consisting of 15 s at 98°C, 20 s at 55°C and 20 s at 72°C and a final elongation at 72°C for 2 min).The purified PCR amplicons were then used as a template for the second PCR (2 min at 98°C, 20 amplification cycles consisting of 15 s at 98°C, 20 s at 67°C and 20 s at 72°C followed by a final elongation at 72°C for 2 min).Primers for the second PCR included unique sequencing barcodes as well as Illumina adapters.Following purification, DNA concentrations were measured and amplification specificity was checked for all samples using gel electrophoresis.
Four negative controls of sterile H 2 O were included during both extraction and amplification that in no case yielded detectable DNA concentrations (based on gel electrophoresis and measured DNA concentrations).The purified PCR products were pooled in equimolar concentrations and sequenced (2 × 150 bp) on the Illumina HiSeq X-ten platform at the Beijing Genomics Institute (BGI, Shenzhen).

| Data processing
Raw sequencing reads were assigned to specific samples (median: 3,345,375 reads/sample), and remaining barcodes were removed.
Due to the fragment size and sequencing technology used, there was no overlap between forward and reverse reads for the 16S fragment (450 bp).For the smaller 18S fragment (260 bp), there was some overlap between forward and reverse reads, but due to deteriorating sequence quality at the reads' ends, we were not able to reliably merge them.Forward reads consistently had a higher sequence quality and therefore were used for all further analyses in QIIME2 (Bolyen et al., 2018).Doing so shortened the overall fragment lengths, which might limit the achievable taxonomic resolution in our analyses.The reads were submitted to three filtering steps using the dada2 denoise-paired plug-in (Callahan et al., 2016).The reads were (i) trimmed and undefined bases removed at a 5 base cutoff and only reads with a minimum length of 150 base pairs (bp) were maintained; (ii) sequencing errors were corrected based on a 2 M read-simulated error model; and (iii) putative chimeras were removed (Table S3).Unless otherwise mentioned, all successive filtering steps and statistical analyses were performed in R v3.5.2 (R Core Team, 2021).Plots were created using the R package "ggplot2" (Wickham, 2009).
A total of 29,365,161 and 36,885,215 reads were assigned to 30,633 and 2748 prokaryotic and eukaryotic amplicon sequence variants (ASVs), respectively (Table S4).Additional filtering included the removal of ASVs with an abundance of more than 20% across the negative control samples (as recommended in Wangensteen et al., 2018) and ASVs that were completely unassigned across the curated 16S and 18S rRNA gene Silva databases (release version 132; Quast et al., 2013) and the NCBI GenBank database release 230 (Johnson et al., 2008).To minimize the effect of PCR bias, the data were further filtered by removing any ASVs that were present in less than two of the four replicates using a similar approach as described in Lange et al. (2015).The four replicates were then combined to one representative sample by adding the remaining read counts of the non-sporadic ASVs to increase the reliability of the diversity metrics (Haegeman et al., 2013).
Taxonomy was assigned by aligning the filtered reads to the two databases mentioned in the previous paragraph.This was done using the plug-in feature-classifier classify-consensusvsearch in QIIME2 (Rognes et al., 2016) and blastn followed by MEGAN (Huson et al., 2016) to find the LCA (last common ancestor).A 97% identity cutoff (of high observed confidence levels in taxonomic classification) was found to be more fitting for these particular datasets, giving accurate and consistent taxonomic classification against the two databases.The taxonomic assignments from the two databases were harmonized to one representative ASV in the abundance tables.After discarding all ASVs that could not be classified at a taxonomy level beyond prokaryotes and eukaryotes, we obtained 11,667 and 1093 ASVs, respectively (Table S5).Rarefaction of all ASVs was performed at a sequencing depth of 2,036,502 and 2,392,669 reads for prokaryotes (including bacteria and archaea) and eukaryotes, respectively (Figure S1).This was done using the rarefy function from the R "vegan" package (Oksanen et al., 2019).Following rarefaction, normalization using the ASV read proportions was used.At this cutoff, full saturation of taxonomic richness was achieved.
To determine community composition across the aquatic environments, relative read abundances were calculated for prokaryotic and eukaryotic phyla (Tables S6 and S7).For illustration purposes and ease of analyses, only phyla present in more than five locations were considered.The barplots corresponding to prokaryotes (i.e., bacteria and archaea) and eukaryotes (Figure 2) were plotted based on the phyla's percent contribution to the total read abundances within each aquatic environment.For two of the most highly abundant eukaryotic groups, Ochrophyta and Arthropoda, taxonomy was further assigned at the family (for Ochrophyta) and class/ subphylum (for Arthropoda) level (Figure S3).

| Data analysis
Alpha diversity (α-diversity) across the aquatic environments was inferred based on two variables: richness (the number of ASVs) and the effective number of species-Shannon-Wiener diversity index (ENS-Shannon; Figure 3).A range of Hill numbers (Chao et al., 2014), including Shannon, Simpson, and Pielou, were calculated using the diversity function from the R package "fossil" (Vavrek, 2011).Nonparametric Spearman's rank correlation coefficients were calculated between the Hill numbers using the cor function from the R package "R stats."Because of their highly correlated values across both prokaryotic and eukaryotic datasets (Figure S2), we chose ENS-Shannon as a representative for all further analyses.The ENS-Shannon index is known to be sensitive to diversity patterns within taxonomic groups and shows no bias due to small sample sizes (Wagner et al., 2018).Variation within each diversity index was calculated using the R function var, and the tests for equality of variances were conducted using the R function levene.testfrom the R package "lawstat."Beta diversity (β-diversity) was calculated as Bray-Curtis dissimilarity values using the function metaMDSdist in the R package "vegan" (Oksanen et al., 2019).Ordinations plots were generated using principal coordinate analysis (PCoA) based on the eigenvectors of the Euclidean distances generated from the reduction of the Bray-Curtis dissimilarity values to principal coordinates.This was done to illustrate the dissimilarity of community composition across all aquatic environments on two-dimensional planes (Figure 4a,b).
To measure the overall average distance of all aquatic environments from the centroids within prokaryotes and eukaryotes, we used the function betadisper in the R package "vegan" (Figure 4c).Statistical difference in means of the normalized distances was calculated Environmental parameters were collected in Nicaragua in 2018 or obtained from previous studies (Table 1; Barluenga & Meyer, 2010;Kautt et al., 2018).We performed Mantel tests based on 9999 permutations of Spearman's correlations to assess the association of β-diversity of prokaryotic and eukaryotic communities with each of the environmental factors (Tables S8 and S9).This was done by generating distance matrices for (i) each environmental factor (excluding geographic distances, see below) based on Euclidean distances using the function dist and (ii) each biological domain (i.e., prokaryotes and eukaryotes) using the function vegdist from the R "vegan" package.
For geographic distance, a matrix based on Haversine distances was generated using the function distm from the R "geosphere" package (Hijmans, 2019).Associations between environmental factors and both richness and ENS-Shannon were estimated by linear regression using the R function lm (Tables S10 and S11), according to the general equation structure: alpha diversity measure (ASV richness, ENS-Shannon) ~ environmental factor (surface area, mean depth, dissolved oxygen, salinity, pH, and age).We performed linear regressions separately since combining multiple explanatory variables and possible interactions in one model was not adequate due to low number of water bodies included in our study.For all statistical analyses, we performed two separate tests: one including all aquatic environments and one for crater lakes only.This was done due to the strong differences between the crater lakes and great lakes (i.e., in α-diversity and in environmental factors such as age) that might have overshadowed the expected subtler differences among crater lakes.
Associations (p < 0.05) between environmental factors and biological communities are illustrated for crater lakes only (Figure 6) while statistics for both subsets are listed in Tables S8-S11.
To obtain information on the extent to which biodiversity is correlated across taxonomic groups, we calculated Spearman's rank correlation coefficients among the subset of phyla that represented more than 1% of the total reads across all aquatic environments (Figure 5).This cutoff eliminated noise from sporadic and less abundant taxonomic groups in the dataset.We tested for differences in correlation coefficients between prokaryotic and eukaryotic phyla using the Wilcoxon rank-sum test (Figure 5b,d;Wilcoxon, 1945) due to non-normal distribution of the data, as determined by a Shapiro-Wilk test (Shapiro & Wilk, 1965).We also tested whether correlation coefficients for the different prokaryotic and eukaryotic phyla were smaller or larger than zero using one-sample Wilcoxon rank-sum tests.

| RE SULTS
To determine the taxonomic composition of prokaryotic and eukaryotic communities for the two great lakes, seven crater lakes and one river across the San Juan drainage basin in Nicaragua, we assigned taxonomy at the phylum level (Figure 2).For prokaryotes, we identified a total of 36 bacterial (mean relative abundance: 98.34% of all sequencing reads across all aquatic environments) and 4 archaeal (mean relative abundance: 1.66%) phyla across all ten aquatic environments (Figure 2a,b, Table S6); only 0.01% of the sequencing reads could not be classified at the phylum level.Proteobacteria (mean relative abundance: 26.38%) and Thaumarchaeota (mean relative abundance: 1.40%) were the most abundant bacterial and archaeal phyla, respectively.For eukaryotes, 28 phyla comprised 91.3% of all sequencing reads (Figure 2a, Table S7), and 8.6% could not be reliably classified at the phylum level.Algae made up the majority of eukaryotic reads (60.88%) with protists, animals, and fungi following at 17.32%, 9.82%, and 3.40%, respectively (Figure 2c).The greatest number of unique ASVs was found for protists (37.97%) followed by algae (24.06%), fungi (14.46%), and animals (10.89%).Among these, Ochrophyta, Ciliophora, and Arthropoda represented the phyla with the highest numbers of unique ASVs, comprising 9.70%, 9.42%, and 6.04% of all ASVs, respectively.In general, the major eukaryotic phyla expected to be found in freshwater environments were well-represented in our dataset (Figure 2c).For the most abundant and sequence-rich algae and animals (Ochrophyta and Arthropoda), we provide a higher taxonomic resolution (Figure S3).The three major arthropod groups (Arachnida, Crustacea, and Insecta) were found across all locations, whereas there was more variation in the geographic distribution of Ochrophyta families (Figure S3).In sum, we successfully classified the majority of sequencing reads to the phylum level, allowing us to describe the geographic distribution of taxonomically variable communities.
We then explored whether prokaryotic and eukaryotic communities differed across the ten aquatic environments using two α-diversity measures: ASV richness and ENS-Shannon, and one β-diversity measure: Bray-Curtis dissimilarity.Prokaryotic and eukaryotic richness and ENS-Shannon were the highest in great lake Nicaragua and Río San Juan and consistently lower in the crater lakes (Figure 3a).Values for great Lake Managua varied for prokaryotes and eukaryotes: Prokaryotic richness and ENS-Shannon were comparable with great Lake Nicaragua, while for eukaryotes, these estimates were as low as those seen in the crater lakes (Figure 3a).
We detected positive correlations between the two α-diversity measures across the aquatic environments for prokaryotic (R = 0.78, p = 0.0075) and eukaryotic (R = 0.7, p = 0.011) communities.
Variances across crater lakes were higher in prokaryotes (richness: We tested for differences in prokaryotic and eukaryotic community composition across the ten aquatic environments based on Bray-Curtis dissimilarity (Figure 3d).Principal coordinate analysis (PCoA) revealed strong differences in the overall community compositions for both prokaryotes and eukaryotes (Figure 4a,b).Across the aquatic environments, distances to the centroid were larger (two-sample t test, p = 0.001) for prokaryotic communities (average Euclidean distance to centroid = 0.612) than for eukaryotic communities (average Euclidean distance to centroid = 0.561; Figure 4c).
Taken together, we detected substantial differences in α-and βdiversity and found that prokaryotic communities were overall more dissimilar across the aquatic environments, which begs the question of what factors are affecting the diversity of these groups.
To test the hypothesis that the variation in prokaryotic and eukaryotic communities across aquatic environments is associated with the variation in abiotic factors, we used linear regression analyses for the two α-diversity measures (richness and ENS-Shannon) and Mantel tests for β-diversity (Figure 6, Tables S8-S11).For prokaryotes, richness (R 2 = 0.831, p = 0.002 and R 2 = 0.860, p < 0.001) and ENS-Shannon (R 2 = 0.78, p = 0.004 and R 2 = 0.880, p < 0.001) increased with age and surface area of water bodies across all aquatic environments (Table S10).When only considering the crater lakes, we found that prokaryotic richness decreased with levels of dissolved oxygen (R 2 = 0.566, p = 0.051) and salinity (R 2 = 0.583, p = 0.046; Figure 6a, Table S11).For eukaryotes, richness (R 2 = 0.86, p = <0.001) and ENS-Shannon (R 2 = 0.88, p < 0.001) increased with surface area across all aquatic environments (Table S10).Analyzing the six crater lakes, eukaryotic ENS-Shannon decreased with the age of crater lakes (R 2 = 0.75, p = 0.027; Figure 6b, Table S11), but this association appears to be mainly driven by the oldest crater Lake Apoyo, which is much older than the other crater lakes.
Differences in prokaryotic community composition (βdiversity) were positively associated with differences in dissolved oxygen levels (Mantel test, R = 0.50, p = 0.011) and salinity (R = 0.56, p = 0.004) across all aquatic environments (Tables S8 and   S9).Considering only the crater lakes, we also detected a positive association between prokaryotic community composition and the difference in dissolved oxygen levels (R = 0.73, p = 0.005) and salinity (R = 0.63, p = 0.006; Figure 6a, Table S9).For eukaryotes, community composition was significantly associated with differences in salinity across all the aquatic environments (R = 0.65, p = 0.005; Table S8) and when considering only crater lakes (R = 0.57, p = 0.005; Figure 6b, Table S9).
We leveraged these differences to test whether prokaryotic and eukaryotic biodiversity was correlated across the aquatic environments.Prokaryotic and eukaryotic richness and ENS-Shannon were not significantly correlated (Spearman correlation, p > 0.05 for both measures; Figure 3b,c).At the phylum level, there was a lot of variation in how the different phyla were correlated (from negatively to no correlation to positively), but overall, mean correlation coefficients across all prokaryotic and eukaryotic phyla were larger than zero for richness and ENS-Shannon (one-sample Wilcoxon rank-sum tests, p < 0.001 for both measures; Figure 5a,c).Further, correlation coefficients for richness ENS-Shannon were higher among prokaryotes compared with eukaryotes (two-sample Wilcoxon ranksum tests, p < 0.001 for both measures; Figure 5b,d); the same result was obtained when only including crater lakes (p < 0.001 for both measures).We tested whether correlation coefficients differed significantly from zero (one-sample Wilcoxon rank-sum tests, p < 0.05) when considering the different phyla separately, that is, by calculating correlation coefficients between a focal phylum and all other phyla.In prokaryotes, correlation coefficients were significantly smaller than zero for no phyla, and larger than zero for 29 and 24 phyla for richness and ENS-Shannon, respectively (out of a total of 40 phyla).In eukaryotes, correlation coefficients were significantly smaller than zero for one and no phyla, and larger than zero for eight and two phyla for richness and ENS-Shannon, respectively (out of a total of 28 phyla).These results confirm that richness and ENS-Shannon indices are overall more positively correlated across phyla in prokaryotes than in eukaryotes.When looking at overall community composition (β-diversity), we found a significant positive association between differences in prokaryotic and eukaryotic communities (Mantel test, R = 0.6822, p = 0.0001; Figure 3d) suggesting that the more different two prokaryotic communities are across two aquatic environments, the more different are the eukaryotic communities as well.

| DISCUSS ION
Understanding which factors contribute to variation in biological diversity of aquatic communities is a common goal of ecology and conservation biology (Reid et al., 2019).Yet, determining aquatic biodiversity is a challenging task (Grey et al., 2018;Gutierrez et al., 2018;Rooney & Azeria, 2015), but eDNA sampling represents an efficient tool to perform large-scale biodiversity monitoring due to its many advantages compared with traditional monitoring techniques (Grey et al., 2018;Lodge et al., 2012;Pawlowski et al., 2020;Thomsen & Willerslev, 2015).While DNA obtained from sediments is now successfully used to reconstruct historical changes in aquatic and terrestrial communities in tropical countries despite numerous challenges (Dommain et al., 2020;Stoof-Leichsenring et al., 2012),

| Biodiversity across heterogeneous aquatic environments
Levels of prokaryotic and eukaryotic α-diversity (richness and Shannon-Wiener diversity index) were consistently higher in the great lakes and Río San Juan compared with the much smaller and younger crater lakes, except for eukaryotic richness in great Lake Managua (Figure 3a).This is congruent with previous studies of teleost fish biodiversity, which showed that the great lakes harbor a large number of species compared to crater lakes (Bussing, 1976;Torres-Dowdall & Meyer, 2021).These results raise the question of why levels of biodiversity differ among these Nicaraguan aquatic environments?The water bodies included in our study vary substantially in an array of variables including age and their physicochemical characteristics (Table 1).Ecological theory suggests that species richness in a given environment is the outcome of origination and extinction of species as well as immigration and emigration rates (MacArthur & Wilson, 1967).These parameters are controlled by different factors, e.g., the geographic area and age of an environment and the extent of connectivity which influences the probability of immigration as well as the availability of ecological niches, as described by the species-area relationship (Preston, 1960).These factors might play a role in the patterns of diversity observed in Nicaraguan aquatic environments.The great and old lakes Nicaragua and Managua have a much larger surface area (8264 and 1024 km 2 , respectively) compared with the small crater lakes (0.2-21.1 km 2 ; Table 1).The great lakes are also much older (500,000 to one million years) than the crater lakes that formed within the last 24,000 years (Bussing, 1976), providing more time and space for colonization events in the great lakes, particularly considering that the crater lakes are completely isolated and have no rivers or streams that flow into or out of them.Further, multiple rivers that differ in water volume, velocity, and productivity flow into the great lakes, increasing the connectivity across water bodies and augmenting the core communities by providing a variety of available habitats.Accordingly, the great lakes harbor higher levels of prokaryotic and eukaryotic biodiversity than the other aquatic environments (Figure 3a), suggesting more productive ecosystems (Gross & Cardinale, 2007;Oehri et al., 2017).
The reduced eukaryotic richness of Lake Managua compared with Lake Nicaragua may be due to the high concentrations of mercury and other toxic substances found in Lake Managua (Lacayo et al., 1991).This confirms previous findings that the diversity of both prokaryotes and eukaryotes within environments is affected by the local conditions, also hinting at the threat of anthropogenic impacts on aquatic biodiversity (Härer et al., 2020;Sun et al., 2019).
Age and surface area of the water bodies were positively associated with prokaryotic α-diversity (richness and ENS-Shannon), while only surface area was positively associated with differences in one measure of eukaryotic α-diversity (ENS-Shannon; Table S10).Similar results have been described for the fish fauna of these lakes: >65 fish species inhabit great lake Nicaragua whereas only between 2 and 19 species inhabit the various crater lakes (Bussing, 1976;Härer et al., 2017;Torres-Dowdall & Meyer, 2021).Even though depth has been shown to be a strong predictor for the diversity of Midas cichlid fishes across lakes (Kautt et al., 2018;Recknagel et al., 2014), this variable does not seem to have a general effect on biodiversity as neither prokaryotic nor eukaryotic diversity was associated with depth.Taken together, age, spatial structuring, and environmental heterogeneity of the great lakes might represent different factors that together can explain the higher levels of prokaryotic and eukaryotic biodiversity, as measured by the two different α-diversity metrics (Figure 3a).Yet, the differences in α-diversity between crater lakes and great lakes might be an underestimation due to our sampling scheme that included one location per water body.Thus, it merely represents a spatial snapshot and additional sampling will most likely detect higher biodiversity levels, particularly in the larger and environmentally heterogeneous great lakes.Hence, we are confident that this major finding of our study (higher levels of α-diversity in great lakes compared to crater lakes) is due to actual biological differences and not due to technical aspects.Additional studies will be necessary to investigate spatiotemporal variation in biodiversity of Nicaraguan lakes, as has been done for other aquatic environments (Lawson Handley et al., 2019).We further speculate that similar analyses on sedimentary DNA retrieved from these environments, historic spatiotemporal inferences can be deduced; this approach has revealed changes in biodiversity over different timescales, ranging from several decades (Ibrahim et al., 2020) to over 100,000 years (Crump et al., 2021).
It should be noted that we were not able to classify particular eukaryotic lineages to higher taxonomic resolution, for example, teleost fishes.One reason might be that the short fragment size we used for our analyses (150 bp) did not contain enough genetic variability to distinguish certain lineages.Recent studies suggest that the volume of water filtered to identify the distribution of tropical fish communities should exceed a certain threshold (>34 L based on recent recommendations by Cantera et al., 2019), which exceeds the volume we used here (i.e., 2 L per water body).However, several studies have successfully detected teleost fishes with volumes of 1-2 L, but these studies were either conducted in temperate regions or in small, artificial environments (e.g., mesocosms, aquaria, and containers), and some used species-specific primers for DNA amplification (reviewed in Rees et al., 2014).All these factors might have increased the potential to capture species diversity in smaller water volumes.In our study, it might be the case that population densities of teleost fish are low in the Nicaraguan lakes, which could make the volume of filtered water a particularly decisive factor for detecting the diversity of fish communities.Furthermore, the resolution provided by the 18S primers used in our study might not be enough to classify chordates to the species level (i.e., different teleost species), but is well-suited for targeting microeukaryotes (Capo et al., 2016).
Therefore, the ability of capturing biodiversity of teleost fishes using eDNA might depend on the primer pair used in the metabarcoding approach, the sequencing technology and the sampling volume.We therefore emphasize the need for additional studies to describe the biodiversity of teleost fishes that should take into account the issues discussed and recommended elsewhere (Rees et al., 2014).

| The effect of abiotic factors on biodiversity
Prokaryotic communities were found to be more dissimilar in their overall community composition across the aquatic environments (β-diversity) in comparison with eukaryotes (Figure 4).The ease of dispersal of prokaryotes facilitates the occupation of ecological niches and colonization of novel environments (Green & Bohannan, 2006;Horner-Devine et al., 2003;Prosser et al., 2007).Thus, we expected prokaryotic communities to be more homogeneous across lakes.Yet, we detected stronger differences in their spatial distribution (Figure 4c), which may reflect differences in local environmental conditions.The Nicaraguan water bodies differ in multiple environmental factors (Table 1) that can alter community composition accordingly (Barnett & Beisner, 2007;Logares et al., 2013;Martiny et al., 2006).The differences in environmental factors can cause habitat diversification by promoting the abundance of habitat specialist taxa (Galand et al., 2009;Pandit et al., 2009).Even though we expected that eukaryotes would show more variation across the aquatic environments due to their limited dispersal across physically detached environments, it appears that prokaryotic communities show comparatively more pronounced differentiation.We hypothesize that this is due to a stronger effect of physicochemical factors on prokaryotic communities than on eukaryotic communities.
Dissolved oxygen levels and salinity affected one α-diversity measure (richness) and β-diversity of prokaryotic communities, and salinity affected β-diversity of eukaryotic communities (Figure 6).
Here, we only focused on the crater lakes as these are more homogeneous in terms of size, age, and geographic isolation, but still differ substantially in their abiotic conditions (e.g., salinity, pH; Table 1).We particularly determined how certain abiotic factors (Table 1), including chemical and physical properties of the water, affect prokaryotic and eukaryotic biodiversity in the seven crater lakes.Salinity can vary greatly across aquatic environments and it is an important determinant of the aquatic biota.Although salinity is often found to have little influence on levels of bacterial species richness (Herlemann et al., 2011;Ji et al., 2019;Wang et al., 2011), it has been shown to strongly affect the overall composition of bacterial communities (Herlemann et al., 2011;Wu et al., 2006;Zhong et al., 2016).Across the Nicaraguan lakes, salinity varies from 0.12 to 4.77 ppt, meaning that all lakes would be considered as freshwater or brackish (for comparison, sea water has around 33 ppt).Thus, our results suggest that the differences in salinity suffice to affect prokaryotic as well as eukaryotic communities.
Congruent with our results, dissolved oxygen levels have been shown to affect bacterial community composition (Zhong et al., 2016) and are negatively associated with bacterial richness (Spietz et al., 2015).While salinity and dissolved oxygen can clearly affect the distribution of biodiversity, we would like to emphasize that other important factors, such as temperature or chlorophyll concentration (Blabolil et al., 2021), were not investigated here and could be included in future studies to obtain a more comprehensive picture of the distribution of biodiversity within and across Nicaraguan aquatic environments.
The fact that we only found significant effects of environmental factors on prokaryotic α-diversity might support the notion that levels of biodiversity in prokaryotic communities across Nicaraguan lakes are largely affected by abiotic, that is, chemical and physical properties of their environments.Taken together, this provides a potential explanation for the stronger differentiation of prokaryotic communities across the aquatic habitats included in this study.

| Covariance in species richness across taxa
The prokaryotic and eukaryotic communities were positively associated in terms of the spatial distribution of community composition (β-diversity; Figure 3d).Yet, we did not detect significant correlations between their overall α-diversity values (Figure 3b,c).
The question of whether biodiversity is correlated among a diverse range of taxonomic groups in certain environments is not only important to better understand the spatial distribution of biodiversity but also for developing conservation strategies and designating protected areas (Flather et al., 1997).Studies that addressed this question obtained mixed results.For example, research on animals and plants from Britain and South Africa revealed that biodiversity hotspots commonly do not overlap for different taxa (Prendergast et al., 1993;van Jaarsveld et al., 1998).In contrast, microbial diversity has been shown to promote diversification through species interaction and niche construction, leading to positive associations in diversity across bacteria (Calcagno et al., 2017;Madi et al., 2020).While prokaryotic and eukaryotic α-diversity was not significantly correlated (Figure 3b,c), we found some interesting patterns when analyzing these datasets separately.α-diversity correlations among phyla were significantly more positive for prokaryotes compared with eukaryotes (Figure 5b,d).A possible explanation for the observed differences is that prokaryotic communities are mainly affected by shared abiotic factors, which might translate to environmentally controlled limits on biodiversity (see previous paragraph for a more thorough discussion on abiotic factors).Since α-diversity was substantially higher in the great lakes and Río San Juan, the observed correlations could be driven by strong differences between great lakes and crater lakes.Yet, when repeating this analysis for crater lakes only, we found the same pattern.It would be of interest to test this higher concordance in the diversity of prokaryotic phyla with a larger number of aquatic environments to determine the generality of these results.

| CON CLUS ION
While determining factors that affect the diversity of aquatic com- (iv) salinity and dissolved oxygen affect both prokaryotic α-and βdiversity, whereas salinity only affects eukaryotic β-diversity across the crater lakes; and (v) there is higher concordance in α-diversity among prokaryotic phyla, suggesting that these lineages are more strongly affected by shared abiotic factors.
Horizontal bar plots representing the percent contribution of each phylum to the total read numbers of prokaryotes (a and b) and eukaryotes (c) within each aquatic environment.Phyla are ordered based on their overall abundance from high to low; only phyla that are present in five or more lakes are shown k i n s o z o a C i l i o p h o r a C e r c o z o a E u g l e n o z o a A p i c o m p l e x a H e t e r o k o n t o p h y t a A m o e b o z o a P e r c o l o z o a B a s i d i o m y c o t a A p h e l i d a B l a s t o c l a d i o t r o t r i c h a C n i d a r i a A n n e l i d a P l a t y h e l m i n t h e s R o t i f e r a C h o r d a t a M o l l u s c a A r t h r o p o d a D i n o a g e l l a t a O c h r o p h y t a C h l o r o p h y t a C r y p t o p h y t a C h a r o p h y t a Animalia Fungi Protista Algae H y d r o g e n e d e n t e s F i b r o b a c t e r e s P A U C 3 4 f L e n t i s p h a e r a e D e p e n d e n t i a e D e i n o c o c c u a -T h e r m u s O m i t r o p h i c a e o t a F u s o b a c t e r i a R o k u b a c t e r i a S p i r o c h a e t e s L a t e s c i b a c t e r i a K i r i t i m a t i e l l a e o t a N i t r o s p i r a e A r m a t i m o n a d e t e s P a t e s c i b a c t e r i a F d o b a c t e r i a V e r r u c o m i c r o b i a B a c t e r o i d e t e s A c t i n o b a c t e r i a P l a n c t o m y c e t e s C y a n o b a c t e r i a P r o t e o b a c t e r i a S y n e r g i s t e t e s C a l d i t r i c h a e o t a E l u s i m i c r o b i a T e n e r i c u t e s C u l i s b a c t e r i a C r e n a r c h a e o t a E u r y a r c h a e o t a T h a u m a r c h a e o t a N a n o a r c h a e o t measured as ASV richness and ENS-Shannon, varies strongly across the ten aquatic environments for prokaryotes and eukaryotes (a).Spearman correlations between prokaryotic and eukaryotic communities for ASV richness (b) and ENS-Shannon (c), and a β-diversity measure (Bray-Curtis distances between all samples) (d).There was only a significant correlation between prokaryotic and eukaryotic communities for β-diversity (d) Principal coordinate analysis plots illustrating the eigenvalues based on Bray-Curtis dissimilarities in unconstrained prokaryotic (a) and eukaryotic (b) community compositions (β-diversity) across the ten different aquatic environments.Across the aquatic environments, distances to the two centroids are larger for prokaryotes (black lines in a) than for eukaryotes (gray lines in b) based on a twosample t test (c), suggesting a higher variation of prokaryotic communities across the aquatic environments, ***p = 0t.testfunction in the R package "stats."The strength and direction of correlations between prokaryotic and eukaryotic αand β-diversity indices were estimated based on (i) Spearman's rank correlations and (ii) Mantel tests using the mantel function from the R package "vegan" (Figure3b-d).

5
Frequency of Spearman's rank correlation coefficients for prokaryotes and eukaryotes based on two α-diversity measures, ASV richness (a and b) and ENS-Shannon (c and d).Mean values for prokaryotes (solid line) and eukaryotes (dotted line) are indicated (a and c).Correlation coefficients are significantly higher in prokaryotes for both measures (b and d).Wilcoxon rank-sum test, ***p < 0.001 Spearman's rho (ASV richness) Association of two α-diversity measures (ASV richness and ENS-Shannon) and one β-diversity measure (pairwise Bray-Curtis distances between all samples) with the environmental factors, shown for prokaryotes (a) and eukaryotes (b) across the crater lakes.We used linear regressions for the α-diversity measures (first two rows in a and b) and Mantel tests for β-diversity (last row in a and b) remained largely unclear how effective eDNA sampling might be in tropical regions to study current patterns of biodiversity across heterogeneous sites.Here, we demonstrate the applicability of eDNA sampling by characterizing and comparing biodiversity across multiple Neotropical aquatic environments that differ in various abiotic factors.This setting proves well-suited to study whether and how the distribution of prokaryotic and eukaryotic aquatic diversity across multiple Neotropical water bodies in Nicaragua is affected by abiotic factors.
munities is important for different biological subdisciplines such as ecology and conservation biology, sampling whole communities has been challenging in the past.Only recently has it become more achievable with emerging technologies like eDNA sampling.By sampling eDNA across multiple Neotropical water bodies in Nicaragua, we collected information on α-& β-diversity of prokaryotic and eukaryotic communities.We successfully identified the taxonomy of prokaryotic and eukaryotic phyla, indicating that eDNA sampling represents an adequate tool to assess prokaryotic and eukaryotic biodiversity in the tropics.Hence, our study provides the first broadscale biodiversity assessment for Nicaraguan lakes using eDNA sampling and highlights its use for capturing biodiversity of tropical aquatic environments.Our main findings constitute that (i) levels of prokaryotic and eukaryotic biodiversity are variable across water bodies and are generally higher in the great lakes; (ii) composition of prokaryotic communities is more dissimilar across lakes than eukaryotic communities; (iii) differences in prokaryotic and eukaryotic community composition (β-diversity) are significantly correlated;