Environmental DNA captures native and non-native fish community variations across the lentic and lotic systems of a megacity.

Globally, urbanization poses a major threat to terrestrial biodiversity, yet its impact on fish diversity is poorly understood, mainly because of surveying difficulties. In this study, environmental DNA metabarcoding was used to survey fish communities at 109 lentic and lotic sites across Beijing, and how environmental variables affect fish biodiversity at fine urban spatial scales was investigated. We identified 52 native and 23 non-native taxa, with lentic and lotic waters harboring both common and habitat-specific species. Water quality strongly affected native fish diversity, especially in lentic systems, but had little influence on non-native diversity. Fish diversity showed little response to urban land cover variation, but the relative sequence abundance of non-natives in lotic waters increased linearly with distance from the city center. Our findings illustrate the complex effects of urbanization on native versus non-native fishes in different aquatic habitats and highlight the distinctive considerations needed to conserve urban aquatic biodiversity.

INTRODUCTION

Urbanization is one of the most substantial anthropogenic processes driving ecosystem and biodiversity changes (, ). Globally, the percentage of people living in cities grew from 3% 200 years ago to 55% by 2018, and it is expected to rise to 68% by 2050 (). City expansions are inevitably accompanied by large-scale land modifications that can lead to extensive and enduring disruption of natural habitats and high rates of native species extinction (). On the other hand, urbanized environments can create habitats favored by introduced non-native species, which may establish and replace the native taxa in urban ecosystems (). Changing biological patterns in response to urbanization, such as reductions in native biota, increases in non-native and invasive species, and biotic homogenization, have been demonstrated in various terrestrial groups, including plants, invertebrates, and vertebrates (, ). Moreover, urbanization often affects biodiversity in a complex and multifaceted manner. For instance, studies of invertebrates, birds, and mammals have demonstrated that, relative to lowly and highly urbanized areas, species richness can peak in intermediately urbanized areas, suggesting that the biodiversity response to urbanization depends strongly on the taxa and spatial scales of investigation (–).

The urbanization process also profoundly alters freshwater habitats and their associated biota. Aquatic habitat changes associated with urban development, including habitat reduction and fragmentation, hydrology modifications, reduced vegetative cover, nutrient and pollutant enrichments, elevated water temperature, and other anthropogenic disturbances, can all alter aquatic community structure (, ). As the dominant vertebrates living in aquatic habitats, fishes are key regulators of nutrient flow and food web dynamics in both freshwater and marine ecosystems (–), and they are often used as biological indicators of water quality and ecosystem integrity (–). In addition, in recent decades, freshwater fishes have experienced rapid population declines and even regional extinctions mainly due to human activities (, ). Therefore, fish communities in urban environments should be monitored because of their pivotal importance in manifesting the effects of anthropogenic pressures on aquatic biodiversity and ecosystems. Understanding urban fish community dynamics can further help to inform ecologically sustainable city planning and management activities.

Studies of fishes in lotic systems (streams and rivers) have typically documented a negative relationship between species richness and urbanization intensity (–). In addition, the urbanization process is often accompanied by an increase in the introduction of non-native fishes for economic, recreational, ornamental, or biological control purposes (, ), and introduced fish species have come to dominate native species in lotic systems in many cosmopolitan areas, such as Paris, France (); Florence, Italy (); and multiple cities in the United States and Brazil (). However, knowledge of the effects of urbanization on fish communities has come almost exclusively from studies of lotic systems, while the effects on fish communities in urban lentic systems (ponds, lakes, and reservoirs) are little understood. Lentic fish species’ taxonomic compositions, ecological adaptations, and response patterns to urban alterations may differ from those of their lotic counterparts, so knowledge of lentic communities is indispensable to fully understanding urban fish biodiversity. Furthermore, in terms of both lotic and lentic habitats, few studies have investigated fish assemblage responses to continuous urban gradients at fine spatial scales, thus precluding a deep understanding of the effects of varying urbanization intensities on this critical group.

The paucity of urban fish data may be largely attributed to the difficulties associated with fish surveys, compared to terrestrial taxa surveys, in urban settings. Traditional fish surveys are predominantly capture-based and invasive (e.g., electrofishing, trapping, and netting) and have limited practicality when it comes to city-based water bodies, which are often small. Recently developed environmental DNA (eDNA)–based analyses use DNA released into the environment (water, soil, etc.) to detect the presence of the originating organisms, thus allowing noninvasive and sensitive biological monitoring (, ). Combined with the use of universal polymerase chain reaction (PCR) primers and next-generation sequencing (i.e., metabarcoding), eDNA analyses facilitate large-scale aquatic biodiversity surveys with unprecedented efficiency and cost-effectiveness (, ).

In this study, we used eDNA metabarcoding to investigate the compositions of fish communities in lentic and lotic water systems across the megacity of Beijing, from its city center to its rural edges, and we analyzed the native and non-native fish diversity response to urban environmental variables in different water types. The city of Beijing was established more than 3000 years ago; it has a long history of urbanization, and it has undergone rapid and extensive expansions and landscape modifications along with a nearly fivefold population increase in the past 70 years. With an area of 16,410 km2 and a population of 21.5 million people (2019 census), Beijing is one of the world’s largest megacities. The city is strewn with man-made, seminatural, and natural lentic and lotic systems (fig. S1). Historical records have documented a rich diversity of native fishes in Beijing (85 species belonging to 12 orders, 21 families, and 65 genera), yet recent traditional surveys of suburban and rural sites indicate a near 50% decline in native species richness and stark contractions of distribution ranges (). No systematic fish surveys have been conducted in the city center area. Using the eDNA approach, we first characterized fish taxon [i.e., molecular operational taxonomic unit (MOTU)] compositions at 109 sites (49 lentic and 60 lotic) across Beijing (Fig. 1 and table S1). We then tested the relationships between native and non-native fish diversity and different environmental factors, particularly water quality indices (WQIs) and landscape variables, in the lentic and lotic systems. Different aspects of fish diversity (qualitative and quantitative α diversity and β diversity) were analyzed to capture the effects of urbanization on various dimensions of biodiversity.

Fig. 1.

Map of Beijing showing the sampling sites (N = 109).

Land cover data were retrieved from 10-m-resolution global land cover maps made in 2017 (). Vegetation includes forests, grasslands, and shrublands.

RESULTS

Fish diversity in Beijing

Following the bioinformatics filtering steps (see Materials and Methods and tables S2 and S3) and standardization of sequencing depth across samples (see Materials and Methods and fig. S2), a total of 75 freshwater fish MOTUs were recovered from 109 sites (Fig. 2, table S4, and fig. S3), with a mean of 18.9 detected per site (SD, 6.6; range, 2 to 37). Of these, 66 (88.0%) were identified to the species level, 5 to the genus level, and 4 to the level of family or above. The MOTUs belonged to nine orders of Teleostei, with Cypriniformes being the most dominant order (55 MOTUs), followed by Gobiiformes (7 MOTUs), and Siluriformes (5 MOTUs) (table S4). Among the 75 freshwater fish MOTUs detected, 52 were native (i.e., documented in historical records of fishes in the Beijing area). The other 23 MOTUs (belonging to five orders, 13 families, and 20 genera) were non-native, including 7 found in the surrounding regions in North and East China, 8 with known distributions in South China, and 8 that do not naturally occur in China (table S4). Four of the foreign taxa [channel catfish (Ictalurus punctatus), North African catfish (Clarias gariepinus), largemouth bass (Micropterus salmoides), and tilapia (Oreochromis spp.)] are recognized as invasive species in China (, ). Non-native taxa were detected at 92.5% of the sites (93.8 and 91.7% of the lentic and lotic sites, respectively).

Fig. 2.

Bar plots showing the numbers of lentic and lotic sites in which each fish taxon was detected.

Native and non-native fishes are each ranked by decreasing number of site occurrences.

Comparison of lentic and lotic fish communities

Overall, 59 (44 native and 15 non-native) and 66 (45 native and 21 non-native) fish MOTUs were found at the lentic and lotic sites, respectively (Fig. 2 and table S4). Of the total 75 MOTUs, 9 (7 native and 2 non-native) were detected only at the lentic sites and 16 (8 native and 8 non-native) were detected only at the lotic sites, suggesting a taxonomic difference in the community compositions of these water habitats. We evaluated the per-site fish taxonomic diversity on the basis of species richness (i.e., the number of MOTUs) and a quantitative index (Shannon diversity), which was calculated using relative reads abundance (RRA) data (see Materials and Methods). Per site, lentic systems supported significantly greater species richness [mean (SD), 21.4 (6.7)] than lotic systems [16.9 (5.9)] (Wilcoxon rank sum test, P < 0.001; Fig. 3 and table S5). However, fish communities from both water types showed similar Shannon diversity values (table S5), indicating little difference in quantitative diversity. Distinguishing native from non-native taxa showed that native species richness was higher per lentic site than per lotic site (P < 0.001), but non-native species richness showed no significant difference between the water types (Fig. 3 and table S5). The β diversity (i.e., areal variation), measured by the average Jaccard dissimilarity (qualitative) between sites, was consistently greater in lotic systems than in lentic systems for natives, non-natives, and total fishes (Fig. 3 and table S5). However, Bray-Curtis dissimilarity (quantitative; calculated on the basis of RRA data) between sites showed more variable patterns, with greater values for non-natives in the lotic systems and greater values for natives and total fishes in the lentic systems (P < 0.001 for all; table S5). Native, non-native, and total fish community compositions all differed between the lentic and lotic sites [permutational multivariate analysis of variance (PERMANOVA), P < 0.05 for all; table S6], but the low R2 values suggested that other environmental and biological factors besides water type also played major roles in structuring fish communities.

Fig. 3.

Comparisons of fish diversity detected at lentic and lotic sites.

Panel A to C show fish α diversity (species richness) and D to F show β diversity (Jaccard dissimilarity) of fish communities. The green and orange dots represent individual data points, the red dots in the center represent the mean, the bars represent the median and the interquartile range, the whiskers represent the 100% percentile range excluding outliers, and the violin plots show the estimated kernel density (frequency distribution). The P values of Wilcoxon rank sum tests are shown.

Native and non-native fish diversity response to environmental variables

We used generalized additive models (GAMs) to investigate the influences of WQIs and landscape variables on the α diversity metrics (species richness and Shannon diversity) of fishes across sites (see Materials and Methods; the summary values of the diversity metrics and environmental variables are shown in Table 1, and the variable selection and model evaluation results are shown in fig. S4 and table S7). While considering all fishes and water types, among the tested WQIs, only chemical oxygen demand (COD) showed a strong negative effect on species richness, whereas total nitrogen (TN), total phosphorus (TP), total dissolved solids (TDS), and dissolved oxygen (DO) each had overall negligible effects (Fig. 4; complete graphs are shown in fig. S5 and statistical results are shown in table S8). Among the landscape variables, only distance from the city center showed a significant impact on fish species richness, while the percentages of vegetation cover and impervious surface, which are the two most common indicators of urbanization intensity, showed little effect (Fig. 4, fig. S5, and table S8). Together, the model explained 55% of the variation in species richness (table S8).

Table 1.

Summary values of diversity metrics and environmental factors of sampling sites.

Distance, distance from the city center; %Veg, percentage of vegetation cover; %IS, percentage of impervious surface; Area, surface area of lentic systems; Width, water width of lotic systems.

	Variable	Mean	SD	Min	Max
Diversity metrics	Species richness	19	7	2	37
	Shannon diversity	4.5	2.4	1.0	14.1
Environmental factors–WQIs	COD (mg/liter)	25	21	0	91
	TN (mg/liter)	3.504	3.430	0.000	14.185
	TDS (mg/liter)	509.0	234.6	168.4	1208.0
	DO (mg/liter)	7.312	4.519	0.009	19.380
	TP (mg/liter)	0.487	0.700	0.000	4.240
Environmental factors–landscape variables	Distance (m)	27,982	23,965	2497	86,432
	%Veg (%)	0.336	0.224	0.002	0.936
	%IS (%)	0.447	0.256	0.012	0.862
	Area (km2)	0.469	0.909	0.003	5.011
	Width (m)	121.75	149.14	6.86	663.68

Fig. 4.

GAMs showing the partial effects of environmental variables on fish species richness.

The tick marks on the x axis represent the data distribution, and the y axis represents the effects of the respective variables. The gray dots represented individual data points, the solid lines are the fitted lines, and the shaded areas represent the 95% confidence intervals. Distance, distance from the city center. ●P < 0.1, *P < 0.05, **P < 0.01, and ***P < 0.001. Results for the less impactful environmental variables, including percentage of vegetation cover (%Veg), percentage of impervious surface (%IS), TN, TP, TDS, DO, lentic system surface area, and lotic system water width, are not shown for presentation clarity. See figs. S5 and S6 and tables S8 to S10 for complete results and GAM splines for Shannon diversity.

Retaining the most influential variables, we applied the GAMs to fish diversity at the lentic and lotic sites separately. The COD remained a strong negative predictor of species richness in both water types, but the effect was more pronounced on lentic than lotic fish species richness, as evidenced by a nearly linear correlation and considerably smaller P values for the former than for the latter (Fig. 4 and table S8). Lotic fish species richness showed a nonmonotonic but overall increasing trend with increasing distance from the city center, whereas lentic species richness did not. Percentage of vegetation cover and TN had little effect on fish species richness in both water types, and neither did surface area in the lentic systems (marginal significance) nor water width in the lotic systems (fig. S5 and table S8).

Next, we analyzed the effects of environmental variables on native and non-native fish diversity separately (Fig. 4, fig. S5, and table S8). Both native and non-native fish species richness varied with COD level and distance from the city center. However, distinct patterns were noted. For instance, native fish species richness decreased with increasing COD across all water types, whereas non-native species richness was only slightly reduced at high COD levels in lentic waters, with no response in lotic waters. In addition, non-native species richness in lotic waters showed a strong increase with increasing distance from the city center (Fig. 4). Compared with the GAM results pertaining to species richness, the Shannon diversity measures typically revealed flattened responses to the environmental variables (fig. S6 and table S9; the statistical significance of each GAM is summarized in table S10).

Relative abundances of native and non-native fishes

Although the GAM results suggested that both native and non-native species richness varied with distance from the city center, there was no significant linear relationship between native or non-native fish diversity (species richness and Shannon diversity) and distance (fig. S7), except for a slight increasing trend of non-native species richness with increasing distance from the city center (R2 = 0.040, P = 0.067). We then investigated the relative contributions of native and non-native taxa to the total communities along the urban gradient by using proportional eDNA sequence abundances (i.e., RRA), which provide a semiquantitative approximation of the relative abundances (biomass or population size) of different taxa within a community (–). We found significantly negative correlations between the RRA of native fishes and distance from the city center across all sites and at the lotic sites but not at the lentic sites (Fig. 5).

Fig. 5.

Correlations of native fish sequence abundance with distance from the city center.

RRA of native fishes at (A) all sites, (B) lentic sites, and (C) lotic sites are shown. The lines indicate fit with a linear model.

Spatial variation of native and non-native fish communities

We further analyzed the correlation between the difference in fish community composition (i.e., β diversity) and the geographic distance between sites. For lentic sites, we found no significant relationships between geographic distance and the qualitative or quantitative β diversity of both natives and non-natives (Table 2 and fig. S8). In contrast, the β diversity of natives, but not non-natives, at the lotic sites exhibited significantly positive relationships with geographic distance (Mantel tests, P < 0.05 for both Jaccard and Bray-Curtis dissimilarities).

Table 2.

Mantel test results of correlations between β diversity of fish communities and geographic distance between sites.

Significant P values are shown in bold.

		Lentic			Lotic
β diversity	Community	R	P		R	P
Jaccard	All fishes	0.069	0.233		0.208	0.008
	Native	0.112	0.114		0.234	0.001
	Non-native	−0.054	0.757		0.060	0.233
Bray-Curtis	All fishes	0.101	0.134		0.200	0.013
	Native	0.103	0.140		0.185	0.022
	Non-native	−0.019	0.584		−0.059	0.876

DISCUSSION

Empowered by eDNA metabarcoding, we expanded the scope of traditional urban fish studies by surveying water bodies of various types and sizes in Beijing, from rural areas to the extremely populated city center. The extensive dataset that we obtained revealed highly resolved fish community compositions at fine spatial scales and thus enabled a robust assessment of fish diversity responses to environmental parameters, thereby providing unprecedented insight into the complexity of the effects of urbanization on aquatic biodiversity. Of all WQIs tested, the COD was the most influential factor on native fish diversity, particularly in lentic waters, whereas impervious surface and vegetation cover, the two most influential urbanization factors shaping the distribution of bird diversity () and land insect diversity (, ), showed little effect on fish diversity. Unexpectedly, the relative sequence abundance of non-native fishes increased along the city-rural gradient in lotic waters. Our results underscore the ideas that the effects of urban modification on biodiversity can be strongly taxa and habitat type specific and that terrestrial, lentic, and lotic biota may show different response patterns within the same urban landscape.

Lentic and lotic systems have distinct hydrological, physiochemical, and ecological features and each provide favorable habitats for various fish taxa, but the former is often neglected in traditional urban fish surveys because of a lack of effective and nondestructive sampling tools. Our eDNA analysis uncovered distinct lentic and lotic fish assemblages, suggesting that surveys focusing on lotic sites alone, as in most conventional urban fish surveys, may underestimate urban fish diversity and that comprehensive assessments of both water types are necessary to capture the full fish diversity of urban landscapes. Per site, the lentic systems had greater α diversity (species richness) than the lotic sites, but between sites, the β diversity (Jaccard dissimilarity) was lower in the lentic systems (Fig. 3), indicating that lotic fish communities were more variable among sites than lentic fish communities. Concordant with these data, the Mantel test results (Table 2) suggested an overall lack of correlation between the β diversity of fishes and geographic distance among the lentic sites.

The hydrological characteristics of different water types inflict different stressors on fish biota. Urban ponds and small lakes often lack efficient renewal, and thus, they are rather static and prone to the accumulation of pollutants, which creates a hostile environment for intolerant native fishes and colonization opportunities for the more tolerant introduced species. On the other hand, streams and rivers are less susceptible to localized pollution discharges because of their flow regime. Corroborating this difference in flow characteristics, the GAM results showed a considerably more prominent effect of water quality (i.e., COD level) on the native fish diversity of the lentic systems compared to that of the lotic systems. The lack of spatial correlation of the β diversity of native fishes in the lentic waters versus the significant spatial correlation in the lotic sites (Table 2) also suggested limited connectivity among urban ponds and lakes. The highly managed nature of city ponds and lakes prevents connection through natural waterways, seasonal flooding, etc., which thus interrupts the metapopulation processes of aquatic biota, which, in the long term, may result in genetic isolation and local extirpation of native fish species. Another disparate result between the GAMs of the lentic and lotic sites regards the relationship between fish diversity and distance: While fish species richness remained relatively constant at various distances in the lentic systems, it showed torturous response curves with an overall increasing trend with farther distances from the city center in the lotic waters (Fig. 4). The heterogeneous ecological conditions and variable anthropogenic interferences at different lotic sites may be an important reason for this result, but it is likely that other yet unidentified factors also contributed to or even determined this pattern. Furthermore, even the full model explained only 55% of the total variance of fish species richness. Many other environmental, anthropogenic, and biological factors can potentially affect the distribution of fish biodiversity. For instance, bank formation, bed substrate, aquatic vegetation, flow velocity, human disturbance, and human-mediated species introduction can also have profound effects on fish assemblages across city water systems. In summary, large gaps still exist in our understanding of fish ecology in urban freshwater ecosystems, thus warranting further investigation.

The high connectivity among lotic water systems facilitates fish dispersal within catchments and across drainages, an effect supported by the significant spatial correlations in the β diversity of native fishes among the lotic sites (Table 2). However, the high connectivity among lotic habitats may also promote the range expansion of non-native species. The significant increase in the relative sequence abundance of non-native fishes in the lotic sites with increasing distance from the city center (Fig. 5) suggests that non-natives may more readily spread and colonize river systems than ponds and lakes. Alternatively, rural rivers may provide favorable habitats for some non-native species (e.g., Rhinogobius brunneus and Zacco acanthogenys) and facilitate their population growth, thus accounting for their high abundances at many lotic sites distant from the city center. Therefore, different water types require specific and tailored management strategies that correspond with particular biodiversity pressures. Native fish conservation calls for strict management of the input of wastewater and other pollutants into city ponds and lakes, and there should be stringent monitoring and control of non-native populations particularly in rivers. Furthermore, it is critical to recognize the importance of connectivity among water bodies and to manage biodiversity in the context of the entire urban water systems ().

Multiple urban stressors act synergistically to exert negative effects on native fish communities, making urban fresh water among the most heavily invaded ecosystems (). Because the traits of non-native fishes are artificially selected to fulfill human needs, introduced species are often characterized by rapid growth, great adaptability, and high tolerance to environmental stress, all of which give them a competitive advantage over their native counterparts (). Studies of freshwater systems with introduced fish species have shown that non-natives have higher tolerance to anthropogenic stress than native fishes and can replace natives in polluted waters (, ). The COD measures the amount of oxygen needed to oxidize the organic matter present in water and is a common indicator of organic contamination by domestic and industrial waste. The WQI strongly inhibited native but not non-native species richness in our analyses, thus suggesting that certain native taxa may be more sensitive to elevated levels of organic input in urban waters and hence be excluded from polluted sites where tolerant species dominate the fish communities. For instance, the minnow Zacco platypus and bitterling Acheilognathus macropterus, native species that are sensitive to wastewater effluents () and hypoxia (), were both found at more than one-third of the sites with low COD values (<20 mg/liter; 58 sites) but at only zero or one site with high COD values (>40 mg/liter; 18 sites) (table S11). In contrast, the three most prevalent non-native species, the gobies Rhinogobius formosanus and R. brunneus and the minnow Sinibrama macrops, all showed similar frequencies of occurrence in high- and low-COD conditions, suggesting that COD had little effect on these species. It should be noted that several of the most common native species in Beijing also showed little response to the COD condition (table S11). The high environmental tolerance and adaptability of these species likely contribute to their prevalence in urban waters, and some of these natives have even become invasive species in other regions [e.g., the stone moroko (Pseudorasbora parva) in Europe and North America ()]. Furthermore, compared with the GAMs of species richness response, the GAMs of Shannon diversity showed flattened responses to COD level. Shannon diversity accounts for the sequence abundance of individual taxa and thus represents a quantitative measure of community diversity. The GAM results suggest that the species sensitive to water quality properties such as COD may be present in low abundance relative to the more tolerant fish taxa. We found no clear fish diversity response to other nutrient and chemical indexes (TN, TP, TDS, and DO). This is not unexpected, as fishes are generally less sensitive to water quality than other aquatic organisms (e.g., macroinvertebrates and diatoms) (, ), and different studies have yielded inconsistent or even conflicting relationships between fish diversity and water nutrient concentrations (). In addition, as Beijing has undergone extensive urban development in the past few decades, the most environmentally sensitive species may have already been eliminated from the city. We caution that different cities may vary considerably in terms of their urbanization histories and land and water modification schemes, which can lead to varying biodiversity responses and should be considered when interpreting observed patterns ().

One of our most unexpected findings was the linear increase of relative sequence abundance of non-natives in lotic fish communities along the city-rural gradient. Previous studies of urban fish assemblages have almost universally shown increased native fish species replacement by introduced species with increasing levels of urbanization (, ). One possible explanation is that, compared to the artificial canals and ditches in city center areas, lotic systems in less urbanized zones have greater habitat diversity and heterogeneity, leading to more suitable environments for introduced species, thereby facilitating their colonization. Non-native species in the study area may have been introduced intentionally or unintentionally via aquaculture transportation, pet trade, religious releases, and water resource engineering projects (, ). For instance, the South-to-North Water Diversion Project, which transports water from the Yangtze River to Beijing and adjacent regions, may have mediated the spread of S. macrops and Odontobutis potamophila (both detected in this study) and potentially other southern China fish species to the receiving water systems. We also noticed that the most common non-native species (e.g., R. brunneus, R. formosanus, and S. macrops; table S4) were all small fishes (typical length < 10 cm), a feature that may provide them greater adaptability and dispersal ability in urban water bodies, which are often small and shallow and frequently constructed with weirs, sluices, and other barriers. Previous studies have shown that certain urban-adaptable fish species (), omnivorous and contaminant-tolerant species in particular, tend to replace specialist and intolerant species (, , ), alter the food web structure, and disrupt ecosystem processes through predation, competition, hybridization, disease transmission, and habitat degradation (). The environmental adaptability and tolerance characteristics of these non-native fishes in Beijing waters await biological and ecological assessments. Further studies of fish physiology, behavioral ecology, life history parameters, adaptive evolution, and population and community biology in urban landscapes will both enhance the mechanistic understanding of fish responses to urbanization and inform aquatic biodiversity conservation practices.

Limitations and methodological considerations relating to eDNA-based surveys

The eDNA method has given rise to new and exciting research avenues. However, eDNA is still a new and developing methodology, and many uncertainties and unknowns remain regarding eDNA ecology and analysis; therefore, caution must be taken during data interpretation (). For one thing, the identification of a fish’s DNA is not necessarily proof of its inhabitance at the sampling site. Apart from shedding by live fish, identification can also result from DNA contamination (e.g., sewage effluents) or assignment errors due to a lack of sufficient interspecific sequence variation. Marine fish eDNA was detected at 10 of the lotic sites in this study, which provides strong evidence of sewage contamination at these sites. However, the small number of lotic sites affected and no such eDNA detected at any of the lentic sites suggested limited influence from such input on the total eDNA data. For the non-native freshwater species detected in our analysis, it is possible that some DNA did not originate from live fishes, although many of these species have also been detected in catch-based surveys (e.g., R. brunneus and O. potamophila) or observed and reported anecdotally (e.g., Acheilognathus imberbis and Oreochromis spp.), supporting the likelihood of their inhabitance in our study area. These newly detected species may represent natural range expansions or human-assisted introductions (see discussion above). In addition, eDNA methods may unveil hidden biodiversity and cryptic species that have eluded traditional surveys (). To cross-examine and evaluate our eDNA data, we compared the species detected in this study with those detected in recent fish surveys using the traditional cage-trapping method in two aquatic systems in Beijing, the lakes in the Old Summer Palace () and the Beiyun River (). The eDNA method and trapping surveys showed largely overlapping fish diversity, but each method also detected unique species (fig. S9). Overall, the eDNA method uncovered a similar or greater number of species than the traditional surveys, but with substantially less sampling time and effort. Our sampling design (three samples per site) was adopted to standardize the sampling effort across sites, but it likely caused the underestimation of species diversity in large water bodies (). With a more thorough and tailored sampling design for specific water systems, we expect that the eDNA-based surveys will demonstrate high efficiency for the detection of fish diversity and that the associated labor and cost will be considerably reduced when using the eDNA method instead of traditional surveys.

The quantitative correlation between the eDNA signal and the actual species abundance is another subject of ongoing debate. Although a number of studies have demonstrated linear relationships between eDNA sequence abundance and fish biomass or population abundance (–), the intrinsic complexities of eDNA ecology warrant a case-by-case validation of this relationship. The most prevalent sequences in our eDNA analysis were consistent with the most common local species [e.g., stone moroko, crucian carp (Carassius auratus), common carp (Cyprinus carpio), and oriental weather loach (Misgurnus anguillicaudatus)], but the quantitative correlation between sequence abundance and fish population size cannot be confidently inferred without formal tests using both eDNA- and catch-based surveys. Despite all these considerations, the detection efficiency of the eDNA assay for individual species was expected to be constant within our study. Therefore, the RRA data were considered comparable across sites to allow for meaningful interpretation of the relative contribution of different fish taxa. However, because the eDNA abundance of individual species was expressed as a percentage ratio at each site, it provided no information regarding the total fish biomass among the sites. Therefore, quantitative measures based on eDNA data may not be directly comparable with those estimated using conventional survey methods and should be interpreted with caution.

Another important technical consideration is the spatial scale of the captured eDNA in relation to its source organism. Once released, eDNA can be dispersed and transported in aquatic environments, meaning that the eDNA detection location may not be closely associated with the actual location of the species. This issue is thought to be particularly problematic in lotic systems because of the water flow carrying eDNA downstream from its site of origin, resulting in an apparent accumulation of biodiversity in downstream eDNA samples (). However, eDNA transport did not likely have a strong effect on the conclusions from our analysis. First, the small water volume and slow flow regime of urban streams and rivers can limit the downstream transport of eDNA. Previous studies have demonstrated that fish eDNA downstream dispersion distances in streams and small rivers are generally within 2 km (–). We placed the sampling sites for the same watercourse at least 3 km apart in longitudinal distance to effectively reduce eDNA carryover between sites. In addition, a high spatial correlation between lotic fish communities is expected when the eDNA transport effect is strong. The lack of spatial correlation between non-native fish assemblages among the lotic sites (Table 2) contradicted this speculation and suggested that eDNA transport did not significantly affect the inferred fish distribution patterns.

Urban ecosystem conservation and eDNA

Cities have created unprecedented ecological problems and habitat damage, but proper management strategies guided by deep understanding of urban biomes can offer remedies to this ecological crisis and can mitigate the conflict between urban living and ecosystem integrity. Shaffer () demonstrated that urban environments can become biodiversity hotspots or even refugia for endangered species. There are few relevant studies relating to aquatic and semiaquatic biodiversity patterns in cities, but they have shown that macroinvertebrates () and amphibians () can exist in high diversity in urban ponds. Integrating biodiversity conservation and ecosystem integrity into urban design and management practices can help create eco-friendly city environments that benefit both ecosystem health and human well-being (, ). Moreover, urban biodiversity provides connections with nature for people living in cities, and it can be a valuable tool to promote nature education and global biological conservation (). The eDNA method enables effective and economical surveys of aquatic biodiversity and also reveals both community dynamics and responses to environmental changes. Some advantages of the eDNA method are its noninvasive nature and the ease of water sampling, both of which make it an ideal choice for public involvement in urban biodiversity monitoring. Citizen science projects using publicly collected eDNA samples have shown promise in surveys of endangered amphibians () and surveillance of invasive species (). With careful study design and method validation, citizen science projects used routinely in eDNA-based fish monitoring efforts may extend our understanding of the consequences of urbanization on aquatic ecosystems and assist in assessing the effectiveness of conservation measures.

In conclusion, fish diversity responses to urban environmental variables differ from those of terrestrial biota, and lentic and lotic systems can exert distinct pressures on native fish communities. Native fish diversity conservation and foreign species control efforts should stress connectivity among water bodies and adopt biodiversity management strategies on the landscape scale for water systems. The noninvasive and sensitive eDNA approach is highly efficacious for examining fish diversity in urban settings and for elucidating the effects of urbanization on fish assemblages at fine taxonomic and spatial resolutions. Researchers can take advantage of eDNA and engage the public to advance our knowledge of urban aquatic biodiversity and safeguard city freshwater ecosystems vital to human life ().

MATERIALS AND METHODS

Water sample collection and environmental factor measurements

Water sampling was carried out in Beijing (39.43°N to 41.05°N, 115.42°E to 117.50°E) (Fig. 1) from April to June 2018 (table S1). Sampling was conducted at 49 lentic sites and 60 lotic sites covering all five river systems flowing through the Beijing area (Juma River, Yongding River, North Canal, Chaobai River, and Ji Canal). A lentic sampling site represented a single lentic water body (surface area, 0.003 to 5.011 km2), from which 1 liter of water samples [a commonly used sample volume for fish eDNA (, )] was collected at three locations spaced roughly equidistant from each other along the bank of the water body. A lotic sampling site represented a segment of a lotic system (surface width, 6.9 to 663.7 m) from which 1 liter of water samples was collected at three locations about 20 m apart from each other. When more than one sampling was conducted for one lotic system, two sampling sites were spaced more than 3 km apart (longitudinal distance). We stored all water samples on ice and filtered them within 8 hours. Measurements of the WQI parameters, including TDS, electrical conductivity (EC), DO, COD, TP, TN, and nitrate-nitrogen (NO3─N), were carried out as described in a previous study ().

The geomorphology of the water systems in the study area is predominately man-made and mostly managed. Water systems in the city center areas, suburbs, and parts of the rural areas are characterized by hardened beds and banks (concrete, slate, or other impervious materials) with thin sediments. Water is generally shallow (<2 m), slow-moving (lotic waters < 1 m/s), and regulated to prevent overflow during the wet season and to maintain minimum water levels during the dry season. The only exceptions are in rural areas where several wetlands and reservoirs have more natural bed substrates and can reach a depth up to 30 m. Therefore, we did not include bed/bank material types, sedimentation, water depth, or flow rate as environmental variables. To assess urban environmental effects on fish diversity distribution, we examined three landscape metrics: (i) percentage of impervious surface (%IS; buildings, roads, and other infrastructure); (ii) percentage of vegetation cover (%Veg; forests, shrubs, and grassland); and (iii) Euclidean distance from the city center, Tiananmen Square (Dist). The first two metrics are commonly used as proxies of urbanization level, and we included the third metric because the development and expansion of Beijing follows a radial pattern characterized by a ring road system with Tiananmen Square at the center, so Dist provided an explicit and intuitive index of urbanization intensity. We obtained Beijing land cover information from high-resolution (10 m) land cover maps () downloaded from http://data.ess.tsinghua.edu.cn/fromglc10_2017v01.html. The percentages of each land cover type (Fig. 1), acquired using ArcGIS 10.3 (), were extracted from a 1-km radius area surrounding each sampling location, and each of those values was calculated as the proportion of the relevant land cover type to total terrestrial area (excluding water areas). The 1-km radius was selected on the basis of previous research on correlations between urban or human impact and biodiversity (, ). Land cover data extracted within a 500-m radius area surrounding each sampling location were also analyzed following the same procedures, and we found highly consistent values for each land cover type with those of the 1-km radius area. Therefore, we only conducted the following analysis using data from the 1-km radius. We averaged the land cover values of the three locations at each sampling site and used those means in subsequent analyses. The Euclidean distances between the city center and the sites were based on the sampling locations’ GPS coordinates calculated with the distm function in the R package geosphere (https://CRAN.R-project.org/package=geosphere) using the Vincenty (ellipsoid) great circle distance method.

eDNA processing and sequencing

Filtration was conducted in a clean laboratory designated for eDNA processing. We filtered each 1-liter water sample and a filtration blank [double-distilled H2O (ddH2O)] through a mixed cellulose esters membrane filter with a pore size of 1.2 μm (Merck Millipore Ltd., Carrigtwohill, Ireland). Filtration units and laboratory tools were treated with a 10% bleach solution (final concentration of hypochlorous acid ~1%) and rinsed with ddH2O twice between samples to remove potential eDNA contamination. Filters were stored at −20°C until DNA extraction.

For DNA extraction, we first finely cut each filter and sonicated it, and then we extracted the eDNA using the DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany) with slight modifications (). A clean filter (an extraction blank) was processed with each extraction batch to monitor contamination.

For the PCR, we used the fish metabarcoding primer set Ac12S (), which targets an ~390–base pair (bp) 12S ribosomal RNA gene fragment and has shown superior and efficient fish diversity detection in Beijing waters (). PCRs of eDNA extracts and negative controls (filtration, extraction, and no-template PCR blanks) were carried out in triplicate, each with uniquely tagged primers to allow for the identification of individual PCR amplicons (). Hence, nine uniquely tagged PCRs were conducted for each sampling site, resulting in a total of 9 × 109 = 981 eDNA PCRs, with an additional 213 negative-control PCRs (filtration, extraction, and PCR blanks). We performed each PCR in a 25-μl total volume that included 5 μl of eDNA extract that was diluted 10-fold (to reduce the PCR inhibitors), 0.2 μM each of the forward and reverse primers, bovine serum albumin (0.4 μg/μl), and 1× Premix Ex Taq (Takara, Kusatsu, Japan), which has ~4.5 times higher fidelity than the standard Taq. The PCR program included an initial denaturation step at 95°C for 10 min; then 45 cycles of 95°C for 30 s, 65°C for 30 s, and 72°C for 1 min; and a final elongation step at 72°C for 10 min (). The PCR products were checked on agarose gels. For the next-generation sequencing, six sequencing libraries were constructed, each containing 132 to 216 PCR products pooled in equal volumes and purified using the EasyPure PCR Purification Kit (TransGen Biotech, Beijing, China). Library construction and sequencing were performed by the Beijing Genomics Institute sequencing service in Wuhan, China, using 2 × 250 bp paired-end sequencing on a HiSeq 2500 System (Illumina Inc., San Diego, CA, USA).

Bioinformatics processing

To process the resulting Illumina sequencing reads, we used programs in the OBITools package (). First, the forward and reverse sequences of the same PCR were aligned using the illuminapairedend program, and sequences with quality scores of <40 were removed using the obigrep program. We used the ngsfilter program to retain sequences with no mismatches in their tags and with a maximum of two mismatches in each primer. Identical sequences were then clustered into unique sequences using the obiuniq program, and sequences with <100 bp or with a total count in the whole dataset of <10 were removed using the obigrep program. We used the obiclean program (abundance ratio = 0.5; sequence difference = 1) to identify and remove putative PCR or sequencing errors (). To minimize cross contamination and tag jumps (), the largest count of a sequence in the negative-control PCRs within a library was subtracted from the sequence count in each eDNA PCR. We used the ecotag program to assign sequences to taxa, and only sequences assigned to the family level or below (genus and species) were retained. We then used the R package lulu () to combine sequences identified as “parent-child” sets with both 97% similarity and nested sample distributions. Information about the sequence read counts after each library processing step is shown in table S2.

We used the Basic Local Alignment Search Tool (BLAST) to search all the unique sequences in the GenBank nucleotide database to make taxonomic (MOTU) assignments. A taxon was assigned to each query sequence on the basis of the Max score and 100% coverage of the sequence. The Max score integrates the rewards for matched nucleotides and penalties for mismatches and gaps and normally gives the same sorting order the traditional e(xpect) value (www.ncbi.nlm.nih.gov/Web/Newsltr/V15N2/BLView.html). For queries that matched only native or only non-native species with Max scores, the species (when only one matched) or the lowest taxonomic level including all matched species (when more than one matched) was assigned. Twenty-two query sequences matched both native and non-native species with Max scores. Because each of those native species also exclusively matched other query sequences with Max scores (i.e., they were independently confirmed MOTUs) while the non-native species did not show independent matches, we assigned the native species to these queries. Sequences with a maximum similarity of <97% (at 100% query coverage) were discarded, and only sequences assigned to fishes were retained for further analysis. Fish MOTUs with read counts <10 or <0.1% of the total reads of each PCR were also removed. Twelve marine fish MOTUs appeared at 10 lotic sites (each MOTU at 1 to 5 sites; table S3), with no particular pattern of geographic distribution. These MOTUs were attributed to DNA contamination from sewage and were disregarded from further analysis. All other MOTUs had known distributions in the study area or possibly occurred there on the basis of their documented distribution ranges, physiological characteristics, or local sighting records. We determined fish species nativeness from historical fish records in Beijing.

Fish community diversity measurements

An average of 53,100 (SD 27,405) sequence reads was obtained from each eDNA PCR following the above bioinformatics processing. To evaluate the sequencing depth effects on the number of detected fish MOTUs, we constructed rarefaction curves for each of the 981 eDNA PCRs with increasing sequencing depth (i.e., number of sequence reads) using the rarecurve function in the R package vegan (fig. S2) (). To standardize the sequencing depth across samples, the sequence reads of each PCR were rarefied to 5000 using the R package GUniFrac (), and 40 of 981 PCRs with <5000 reads were discarded. This sequencing depth was sufficient to capture the detected fish MOTUs across samples (fig. S2).

Fish community data were analyzed using both incidence-based and abundance-based metrics across sampling sites. The incidence-based matrix was summarized as the presence (i.e., detected in one or more PCRs of each site) or absence (i.e., not detected in any of the nine PCRs of each site) of each MOTU, and the abundance-based matrix was summarized as the RRA (i.e., percentage of reads for each MOTU, averaged across the nine PCRs of each site; table S4).

We estimated taxonomic diversity (i.e., α diversity), the classical measure of biodiversity, using species richness (i.e., the number of fish MOTUs; equals to the Hill number when q = 0) and Shannon diversity, calculated as the exponential of classic Shannon entropy [exp(−∑pilogpi), where pi is the RRA of each MOTUi at each site; equals to the Hill number with q = 1) with the R package hilldiv (). The Hill numbers provide a unified statistical framework for easily interpretable and comparable diversity measures and have been increasingly adopted in DNA metabarcoding analyses (, ).

We measured the differences in fish community compositions among different sampling sites (i.e., β diversity) both qualitatively and quantitatively using the presence/absence-based Jaccard coefficient and the sequence abundance (i.e., RRA)–based Bray-Curtis coefficients, respectively, calculated using the vegdist commands in vegan. Both the Jaccard and Bray-Curtis coefficients range from 0 to 1, with a value of 0 indicating completely identical community compositions and a value of 1 indicating communities without any shared common taxa.

Statistical analyses

We used Wilcoxon rank sum tests to examine statistical differences in the diversity metrics (i.e., α diversity and β diversity) between the lentic and lotic sites. To test for differences in fish community compositions between the lentic and lotic sites, we performed PERMANOVA using the Jaccard and Bray-Curtis coefficients as the response variables, respectively. The PERMANOVA was conducted using the adonis function in vegan with 999 permutations.

As biodiversity responses to the environment can be nonlinear and nonmonotonic, we used GAMs () to analyze the relationships between selected environmental factors and fish diversity metrics. GAMs are highly flexible models in which the effects of the explanatory variables on the response are defined by data-driven smooth functions rather than a priori (e.g., linear) patterns. We implemented the GAMs by using the function gam with default settings in the R package mgcv (). To construct the GAM for the diversity metrics, we first used the pairwise Spearman’s rank correlation to assess the collinearity between the predictor variables (seven WQIs and three landscape factors; fig. S4), and then we removed one of each pair of highly correlated variables (correlation coefficient |r| > 0.7). This resulted in EC (correlating with TDS) and NO3─N (correlating with TN) being removed. Although Dist and %IS were highly correlated (r = −0.79), they may each represent important and different urbanization characteristics, so we included both in the full model. We ran the full GAM model, which included the three landscape variables (Dist, %Veg, and %IS) and five WQIs (COD, TN, TP, TDS, and DO) as well as additional GAMs using 25 combinations of one to seven of these variables, and then calculated the Akaike information criterion (AIC) values of the models using the function AIC in the R package stats (). The full model was among the best models of all the tested GAMs, according to the AIC values (ΔAIC < 2 compared with the lowest value) and the percentages of deviance explained (table S7). Therefore, we applied the full model to fish diversity (natives and non-natives separately and combined) in all water types. Next, we constructed separate GAMs for lentic and lotic waters. To avoid overparameterizing, we retained only the variables that showed apparent effects (P < 0.1) for all water types in the full model (Dist, %Veg, COD, and TN) and added an additional variable of water body size (surface area for lentic sites and surface width for lotic sites). A Gaussian distribution with an identity link function was applied, and thin-plate regression splines were chosen as smoothers for the separate variables. We evaluated model fit by examining the percentage of deviance explained and residual distribution. Linear relationships between distance and the diversity metrics of native and non-native fishes and between distance and the RRA of native fishes were analyzed with a linear model generated by using the lm function with default settings in the R package stats ().

In addition, we conducted Mantel tests in vegan using the mantel function with 1000 permutations to assess the correlations between fish community dissimilarities (β diversity) and the geographic (Euclidean) distances among sampling sites in lentic and lotic systems. All statistical analyses were conducted in R version 4.0.3 ().