Successive accumulation of biotic assemblages at a fine spatial scale along glacier-fed waters

Summary Glacier-fed waters create strong environmental filtering for biota, whereby different organisms may assume distinct distribution patterns. By using environmental DNA-based metabarcoding, we investigated the multi-group biodiversity distribution patterns of the Parlung No. 4 Glacier, on the Tibetan Plateau. Altogether, 642 taxa were identified from the meltwater stream and the downstream Ranwu Lake, including 125 cyanobacteria, 316 diatom, 183 invertebrate, and 18 vertebrate taxa. As the distance increased from the glacier terminus, community complexity increased via sequential occurrences of cyanobacteria, diatoms, invertebrates, and vertebrates, as well as increasing taxa numbers. The stream and lake showed different community compositions and distinct taxa. Furthermore, the correlations with environmental factors and community assembly mechanisms showed group- and habitat-specific patterns. Our results reveal the rapid spatial succession and increasing community complexity along glacial flowpaths and highlight the varying adaptivity of different organisms, while also providing insight into the ecosystem responses to global change.


INTRODUCTION
Glaciers hold the main freshwater resources on Earth and play pivotal roles in regulating regional hydrology and climate. 1,27][8] Therefore, there is an urgent need to understand the ecological processes of these systems to predict the impacts of glacial retreat on diverse biota and ecosystem dynamics.
Glacier-fed waters are extreme, dynamic, and heterogeneous systems.From an ecological perspective, the low temperature, low atmospheric oxygen, high ultraviolet radiation, frequent disturbance, and ultraoligotrophic status of glacier-fed aquatic systems stress organisms in many ways and create strong environmental filtering for the taxonomic lineages and biological traits of those living in such habitats. 6,9hose environmental stressors are particularly prominent near the glacier terminus and result in biodiversity ''coldspots'' that have low species numbers and productivity. 8,10As glacial influences weaken with distance along the flowpath, glacier-fed waters progress through diverse hydrological and physiochemical features that allow them to ultimately support a highly diverse community of organisms. 6,11Indeed, different biota form spatially distinct communities in and around glacial waterways, thus demonstrating how the variable inhabitability of glacier-influenced environments correlates with the distinct adaptivity of specific taxa. 12While the ways different biological groups assemble and organize along glacier-fed waters are not well understood, the sensitivity and tolerance of organisms to abiotic environments and their distinct requirements for biotic conditions are known to vary. 13,14Resolving the biodiversity and community compositions along glacial waterways will provide vital information about the processes and drivers of community structures in glacier-influenced aquatic ecosystems.Such knowledge can also aid understanding of biological community organization in other extreme environments 15 and help elucidate taxa-specific responses and ecosystem alterations under climate change. 6he traditional methods of performing comprehensive biodiversity assessments are daunting tasks fraught with logistical obstacles because of the need for diverse survey methods and extensively varied expertise in morphotaxonomic identification across vastly different biological groups.The already strenuous task of traditional biodiversity surveys is made even more difficult when coupled with the remote and harsh environmental conditions of most glaciers.Therefore, most studies of glacier-influenced biota have addressed a single biological group (e.g., bacteria, arthropods, and fish).However, recent technological advances in environmental DNA (eDNA) methodology have allowed non-invasive, robust, and economic surveys of diverse organisms via the detection of trace amounts of their DNA in environmental

DNA sequencing and taxonomic assignments
We used four metabarcoding primer sets to amplify eDNA extracts, with each targeting a specific biological group's gene regions: the cyanobacterial 16S rRNA gene (CYA), 25 the diatom rbcL gene (708F), 26 the invertebrate COI gene (BF), 27 and the vertebrate 12S rRNA gene (Tele02) 28 (Table S2).We obtained 127,474,482 raw reads from all sequencing libraries, of which 6,955,873 (cyanobacteria 737,282; diatoms 4,482,675; invertebrates 35,206; vertebrates 1,700,710) remained after stringent quality filtering and bioinformatics processing.The sequences were clustered into 642 operational taxonomic units (OTUs); 125 OTUs were assigned to cyanobacteria, 316 OTUs to diatoms, 183 OTUs to invertebrates, and 18 OTUs to vertebrates (Table S3).The rarefaction curves demonstrated that the OTU accumulations in most of the samples reached an asymptote (Figure S1), indicating that those OTUs sufficiently represented their respective assemblages.Subsequently, the resulting OTU reads were normalized to generate quantitative (i.e., relative read abundance [RRA]) data for further analyses (Data S1).
Among the 125 cyanobacterial OTUs, 90 (72.0%OTUs, 82.1% RRA) were classified at the family or lower taxonomic level (Table S3), and 16 families within 12 orders were identified (Table S4).Of the 316 diatom OTUs, 245 (77.5% OTUs, 88.2% RRA) were successfully classified to the family or lower taxonomic level, and those taxonomic assignments revealed 20 families belonging to 10 orders.Of the 183 invertebrate OTUs, 135 (73.7% OTUs, 63.4% RRA) were classified at the order or lower taxonomic level, and a total of 21 orders within 14 classes were identified.For vertebrates, all 18 OTUs were successfully assigned to the family or lower taxonomic level, and 11 families within eight orders were identified.

Multi-group community compositions at proglacial and lake sites
Cyanobacteria were detected at all sampling sites, with an average 30.5 (range 10-57) OTUs recovered per site (Figure 2; Table S5).At the family level, the cyanobacterial assemblages were dominated by Cyanobiaceae (mean RRA: 54.7%), followed by Leptolyngbyaceae (15.8%) and Pseudanabaenaceae (8.3%), across the sampling sites (Figure 3).Notably, 17.9% of all detected cyanobacterial sequences could not be assigned to a specific cyanobacterial family and were categorized as ''Unassigned.''There were noticeable differences in the relative contributions of various cyanobacterial families between the proglacial and lake samples, with Leptolyngbyaceae and Pseudanabaenaceae predominating in the proglacial samples, whereas Cyanobiaceae was most prevalent among family-assigned OTUs in the lake samples (Figure 3).
Diatoms were recovered at all lake sites and in the glacier-fed stream (PL03-PL05), but not at the first two proglacial sites near the glacier terminus (PL01 and PL02).While the taxonomic richness of diatoms was also low in the glacial stream (4, 3, and 26 OTUs detected at PL03, PL04, and PL05, respectively), it was greater in each lake sample (OTU range 42-136, mean = 95.9)(Figure 2).Across all sites, the diatom assemblages were dominated by Fragilariaceae (mean RRA: 49.4%), followed by Stephanodiscaceae (16.3%) and Cymbellaceae (13.7%).Proglacial and lake samples again showed different compositions of diatom taxa, and OTUs unassigned at the family level yet again accounted for large proportions of the proglacial assemblages (Figure 3).Among the family-assigned OTUs, Cymbellaceae dominated at the first proglacial site where diatoms were detected (PL03), while the proportion of Fragilariaceae gradually increased toward the downstream sites.In contrast, Fragilariaceae was the overall most abundant family across the lake samples, but Stephanodiscaceae increasingly dominated the lower lake sites close to the outflow.
Invertebrates were detected in all lake samples, but they were detected only at the first and last proglacial sites (PL01 and PL05).As with diatoms, the invertebrate taxonomic richness was low at the proglacial sites (6 and 3 OTUs found at PL01 and PL05, respectively), but it was generally higher in the lake samples (OTU range 2-62, mean = 27.2) (Figure 2).The assigned invertebrate taxa at the order level were dominated by Lepidoptera (mean RRA: 32.8%), Diptera (16.8%), and Tubificida (5.7%) across all sites (Figure 3).The only OTUs that were assigned to the order or lower levels at the proglacial sites were Odonata at PL01 and Diptera (family Chironomidae) at PL05.However, OTUs unassigned at the order level accounted for an average 36.6% of invertebrate sequences across sites and were prevalent in both proglacial and lake assemblages.
No vertebrate sequences were recovered from any of the proglacial samples, but they were found at all lake sites (3-8 OTUs per site; Figure 2).The dominant taxa within the vertebrate assemblages were fish (6 OTUs), specifically the carp family Cyprinidae (subfamily Schizothoracinae in particular; mean RRA = 55.1%) and the stone loach family Nemacheilidae (genus Triplophysa; 15.0%), which are prevalent in the region's fish assemblages.The six mammalian OTUs consisted mostly of livestock (cattle Bos taurus, yak B. grunniens, horse Equus caballus, and sheep Ovis aries), and the six detected avian OTUs included both locally occurring wild birds (blood pheasant Ithaginis cruentus, redbilled chough Pyrrhocorax pyrrhocorax, white-breasted waterhen Amaurornis phoenicurus, pigeon Columbia spp., and thrush Turdus spp.) and a domesticated species, the chicken Gallus gallus (Data S1).These findings demonstrate the ability of eDNA analysis to capture both aquatic and terrestrial species' signals in water samples that may be accumulated via direct deposit and surface runoff from surrounding areas.Domesticated species accounted for large proportions of vertebrate OTUs at all lake sites (Figure 2), revealing the considerable influence of human activities around the lake.
Pearson correlation analyses revealed that the a diversity (i.e., the number of OTUs detected at each site; Table S5) of all biological groups increased significantly as the collection sites progressed from the proglacial headwater to the lower lake (Pearson: r = 0.64, p = 0.008 for cyanobacteria; r = 0.63, p = 0.036 for diatoms; r = 0.70, p = 0.026 for invertebrates; r = 0.58, p = 0.019 for vertebrates).
The degree of water eutrophication, as indicated by the diatom-based indicator Pampean Diatom Index (IDP) (see method details), also exhibited a similar increasing trend with the direction of flow (Pearson r = 0.90, p < 0.001).The low IDPs (0.5-1.1) of the proglacial sites where diatoms were detected indicated good water quality with low nutrient levels, but the higher IDP values (i.e., >2) of the lake samples, particularly the lower lake sites (RW08-RW11) (Figure 2), indicated poorer water quality with higher nutrient levels. 29The higher IDP values at the lake sites compared to the proglacial sites were also correlated with strong eDNA signals of domesticated animals in the lake samples.Only one lake site (RW07), located near the inflow of another glacial stream (Figure 1), had a low IDP.

Differences between proglacial and lake communities
The non-metric multidimensional scaling (NMDS) analysis of the cyanobacteria, diatoms, and invertebrates revealed a consistent clustering pattern wherein the assemblages from the proglacial and lake sites were clearly separated (Figure 4).Furthermore, the cyanobacteria The taxonomic levels displayed are based on more than 70% of the OTUs identified for each group having taxonomic information at that level.The figures depict the 10 most abundant families of cyanobacteria, diatoms, and vertebrates, and the 10 most abundant orders of invertebrates.Low abundance families/orders were collapsed into ''Others.''The OTUs that could not be assigned to a single family (for cyanobacteria and diatoms) or order (for invertebrates) were grouped and are shown as ''Unassigned.''PCR results for certain proglacial sites were excluded due to low read abundance after sequence filtering, thus explaining the blank sections on the plots.
We more deeply investigated the differential representation of individual OTUs between the proglacial and lake assemblages to identify species that showed strong habitat specificity and were thus main drivers of community turnover.For cyanobacteria, Cyanobacteriia_OTU3, Leptolyngbyaceae_OTU6, and Cyanobacteriia_OTU5 were more abundant in proglacial samples in comparison with lake samples, while Cyanobium_PCC-6307 showed the opposite trend (Figure 5).Further analysis based on the similarity percentage (SIMPER) showed that Cy-anobium_PCC-6307 accounted for most of the dissimilarity between the proglacial and lake assemblages (contribution = 43.8%;Table S6).
For diatoms, Cymbellales_OTU11 and Encyonema minutum were more abundant in the proglacial samples, whereas the Fragilaria acus/ radians complex, Stephanodiscus hantzschii, and Cymbella neocistula were more abundant in the lake samples (Figure 5).SIMPER results showed that Cymbellales_OTU11 and the Fragilaria acus/radians complex accounted for most of the dissimilarity between proglacial and lake assemblages (Table S6).
For invertebrates, Insecta_OTU47, Chironomidae_OTU8, Insecta_OTU88, and Insecta_OTU57 were the top four most abundant OTUs in proglacial samples.Lepidoptera_OTU2 and Insecta_OTU4 were found only in lake samples and were the two most abundant lake taxa (Figure 5).Lepidoptera_OTU2, Insecta_OTU47, and Chironomidae_OTU8 were the top-ranked OTUs responsible for dissimilarity between the proglacial and lake assemblages.
The compositions of different biological groups showed distinct correlational patterns with environmental variables.Longitude, latitude, and ammonium nitrogen (NH 4 + -N) were the main environmental factors correlated with cyanobacterial assemblage composition (Mantel's r = 0.58-0.94,all p < 0.05).Longitude (r = 0.72, p = 0.025) and total carbon (r = 0.58, p = 0.008) were significantly correlated with diatom assemblages.Invertebrate assemblages were significantly associated with six of the 12 environmental factors, including longitude, latitude, dissolved organic carbon, dissolved organic nitrogen, NH 4 + -N, and total phosphorus (r = 0.68-0.95,all p < 0.05).None of the evaluated environmental factors had significant relationships with vertebrate assemblages (Figure 6).

Ecological assembly processes of different groups
The community assembly mechanisms of different biological groups were assessed using iCAMP, which revealed that the relative contributions of the community assembly processes (selection, dispersal, and drift) varied among different habitats and groups (Figure 7).For cyanobacteria, community assembly at proglacial sites was predominantly influenced by drift (69.7%), followed by homogeneous selection (25.2%), whereas homogeneous selection (66.6%) had a larger impact than drift (30.7%) at lake sites.For diatoms at lake sites, the community assembly was driven primarily by drift (58.9%), followed by dispersal limitation (24.7%).The community assembly of invertebrates at lake sites was influenced primarily by homogeneous selection (50.6%), followed by drift (26.5%) and heterogeneous selection (14.5%).

Biodiversity composition and community build-up in glacier-fed waters
This study is the first, to our knowledge, to report the sequential accumulation of various biological groups as the distance from a glacier increases, together with a gradual increase in community composition complexity along glacial meltwater flowpaths, within a small geographic range ($1.8 km from PL01 to PL05).Our findings highlight the spatial heterogeneity of biodiversity in glacier-fed waters, while illuminating the fine-scaled ecological succession of biological communities under varying levels of glacial influence.Furthermore, our results reveal the differences in environmental tolerance and adaptability among biological groups and taxa.
Cyanobacteria were the only group detected at all sites across the sampled proglacial waters.The pervasiveness of cyanobacteria was expected, because they are known to constitute an important component of microbial communities in cryosphere habitats, such as glaciers and Arctic soils. 30,312][33][34] We found that the Leptolyngbyaceae and Pseudanabaenaceae families were most prevalent in proglacial samples, which contrasted with the Cyanobaceae-dominated lake assemblages (Figure 3).The remarkable turnover in assemblage composition suggests that divergent environmental filtering effects are exerted on cyanobacteria in different glacier-influenced habitats, while revealing the ecological adaptivity of each habitat's characteristic taxa.Leptolyngbyaceae was the most prevalent cyanobacteria in supraglacial sediments/soil from the Arctic archipelago of Svalbard 35,36 and in glacier cryoconite from the Tien Shan Mountains in China. 37Pseudanabaenaceae was one of the most common cyanobacteria families in cryoconite from polar and mountain glaciers around the globe 38 and in soils from glacial forefields in both Alaska and Peru. 39Our finding that Leptolyngbyaceae and Pseudanabaenaceae were abundant in glacial meltwaters on the Tibetan Plateau suggests universally convergent selection of cyanobacteria lineages in glacial ecosystems.Because glacier meltwater may harbor autochthonous bacteria, as well as allochthonous microbial communities that originate from the surrounding snow, ice, rocks, and soil, 40 further investigations of microbial community dynamics could shed light on potential ecological linkages between different glacial habitat compartments.
Diatoms, the second most omnipresent group in this study, were found in all sampled waters except for two proglacial sites near PL's terminus.Few studies have investigated diatoms in meltwaters at the glacier terminus.While diverse diatom taxa have been found in glacier cryoconite holes, their low abundance in that habitat is suggestive of allochthonous deposition rather than habitation by biological The inner bar graphs represent the OTU relative read abundance (RRA) in the proglacial (blue) and lake (green) samples.The outer ring indicates the mean OTU RRA across all sites, and the names of the OTUs that contributed R5% to differences between the glacier and lake assemblages are shown in the color corresponding to their more prevalent habitat.communities. 41,42Diatom taxa diversity gradually increased in the meltwater stream with distance from the glacier and reached high richness downstream in Ranwu Lake, suggesting that the harsh physiochemical properties and ultraoligotrophic conditions at PL's terminus may inhibit diatoms' stable inhabitance.Interestingly, several dominant diatom taxa in our stream samples (i.e., Encyonema minutum, Hannaea arcus, and the genera Cymbella and Fragilaria; Figure 5) were also prevalent in Canadian glacier-fed streams, 43,44 suggesting that glacial diatoms undergo conserved environmental filtering around the globe.Diatoms often dominate freshwater algal communities, and they are the main eukaryotic primary producers in alpine fluvial ecosystems. 45,46Hence, high levels of diatom richness and abundance may be prerequisites that allow organisms occupying higher trophic levels to occur in nutrient-poor glacial waters.
Aqueous eDNA can originate from both water and land-dwelling biota.The dominant invertebrate taxa recovered at the two proglacial sites, Odonata (dragonflies and damselflies) at PL01 and Chironomidae (non-biting midges) at PL05, both have aquatic (larvae) and terrestrial (adult) life cycle stages; therefore, the origin of the eDNA for those taxa could not be determined.The proglacial site PL01 was determined to be an unlikely habitat for carnivorous Odonata larvae because it had a particularly low water temperature and harbored few aquatic animals.Therefore, the Odonata eDNA at PL01 was most likely deposited by flying adult individuals.The lack of invertebrate signals at most proglacial sites also supports their unsuitability as invertebrate habitats.However, Chironomidae detected at the downstream site PL05 may represent stream-dwelling larvae, as morphology-based studies have often found various Chironomid larvae dominating benthic macroinvertebrate assemblages in glacial rivers across Europe. 14,47,48o vertebrates (fishes, birds, and mammals) were detected at proglacial sites, but they were detected at all lake sites, confirming that glacial headwaters, with their poor nutrient levels and low faunal abundance, are unfavorable habitats for large-bodied species.Overall, the sequential addition of different biological groups and increasingly complex community organization along glacier-fed waters reflect the environmental tolerances of different taxa and the ecological drivers of biodiversity that likely encompass both abiotic filtering and biotic interactions. 49

Environmental drivers and community assembly mechanisms
The assemblage composition of different biological groups showed distinct environmental correlation patterns across sites (Figure 6).Cyanobacteria, diatoms, and invertebrates responded to both spatial and water nutritional factors, but the specific patterns varied.The assemblage compositions of cyanobacteria, diatoms, and invertebrates were all correlated with longitude.As environmental variables were predominantly measured at lake sampling sites, and the overall water flow direction in Ranwu Lake is longitudinal, this correlation indicated that spatial distribution variations in those groups may be associated with habitat characteristics, such as hydrological and substrate conditions, at different localities within the lake.Cyanobacteria and invertebrate assemblages were additionally correlated with latitude, suggesting that their assemblages may be more susceptible to the influences of subtle habitat changes in comparison with diatom assemblages.Different groups also showed divergent correlations with water nutrients.While cyanobacteria and diatom compositions were related to a single nitrogen metric and a single carbon metric, respectively, invertebrate assemblages were found to be associated with four metrics of carbon, nitrogen, and phosphate.1][52] In contrast to the microbial and invertebrate communities, vertebrate assemblage composition showed no correlation with spatial or physiochemical water factors, likely because many of the detected vertebrate species were terrestrial animals with strong locomotion ability, which were less affected by spatial factors or water properties.Similar community response patterns among biological groups have been observed for other glacier-fed running waters, where water eDNA-detected microbial and invertebrate communities showed substantially stronger correlations with environmental (spatial, hydrological, and water physiochemical) variables in comparison with vertebrates. 12isentangling the relative contributions of various ecological processes to community assembly is crucial for understanding drivers of biodiversity distribution patterns and forecasting future ecosystem trends in response to environmental changes. 53,54We found disparate mechanisms for cyanobacterial community assembly in different habitats (proglacier vs. Ranwu Lake) and for different biological groups in the same habitat (Figure 7).As the responses of different cyanobacteria taxa to environmental factors can vary greatly, and the proglacial and lake sites showed significant differences in cyanobacterial assemblage composition and environmental conditions, it is not surprising that the community assembly mechanisms for cyanobacteria differed between the two habitats.Cyanobacteria from Ranwu Lake showed primarily homogenous selection-driven assembly, perhaps due to the harsh, yet relatively uniform, environmental conditions in the lake, which may have led to strong and similar forms of selection pressure across the lake.Similarly, homogeneous selection has been shown to dominate microbial community assembly in other extreme and energy-restricted systems, such as glacier-fed streams, 55 an oligotrophic gyre, 56 and deep-water sediments. 57However, cyanobacterial assemblages at the proglacial sampling sites appeared to be an exception to this general pattern.The proglacial sites included a number of localities that showed high habitat heterogeneity, in which water hydrology and physiochemical properties changed rapidly as the distance from the glacier terminus increased.These changes may have led to considerable variations in the cyanobacterial assemblage compositions across the proglacial sites (see also Figure 4) and could have confounded the assembly mechanism analysis.Thus, the contributions of selection and dispersal were both low, while those of drift and other neutral processes dominated cyanobacterial assembly at proglacial sites.
In comparison with prokaryotic biota, the community assembly of eukaryotic biota has been less studied.Limited evidence has shown that, consistent with our result, ecological drift (neutral processes) primarily accounted for the community assembly of coastal phytoplankton 58 and diatoms in a glacier-fed stream. 12Community assembly of diatoms in Ranwu Lake was mainly explained by drift, suggesting that the influences of environmental selection and dispersal on diatoms were limited.Of all biological groups assessed in our study, diatoms were the group with the highest number of OTUs detected.Different diatom species may have distinct environmental responses and dispersal ability, which lead to different distribution patterns within the lake.Therefore, the importance of deterministic processes (selection) was relatively low for diatoms in the cross-lake analysis, while community assembly appeared to be governed by stochastic processes (drift and dispersal).Few studies have examined invertebrate community assembly mechanisms in glacial habitats.Our analysis showed that invertebrates in the lake, similar to cyanobacteria from the same sites, had homogeneous selection-dominated assembly, which was indicative of a strong effect of environmental filtering on this group.Further investigations of other systems are needed to determine whether our findings represent a common pattern for invertebrates in high-altitude lakes.As the climate continues to warm and glaciers melt at a correspondingly faster rate, the relatively influences of different factors on the species composition and community assembly mechanisms of various biota inhabiting glacial environments will likely change. 59More in-depth study is necessary to fully understand the pattern and impact of these alterations.

Knowledge gaps and future studies
Compared with those recovered from the lake sites, substantially larger proportions of OTUs detected in the proglacial samples could only be assigned to cyanobacteria, diatoms, and invertebrates at higher taxonomic levels, while they matched no existing reference sequences with high identity (Figure 2), indicating that the proglacial biodiversity on the Tibetan Plateau is considerably underrepresented in the available sequence databases.In addition to the high endemism characteristic of Tibetan Plateau biodiversity, 60,61 the uniqueness of the glacier biota has been shown to be particularly high at the terminus zone. 10The extreme environmental conditions of the glacial terminus zone and the low habitat connectedness between glaciers likely give rise to specialized biota that are unique to each glacier, and challenges associated with assessing glacial fronts and obtaining biological samples further complicate the difficult task of constructing comprehensive glacial biodiversity databases.For instance, by sequencing metagenomes and cultured isolates from 21 Tibetan glaciers, Liu et al., 2022 62 found that 88.3%-100% of 968 detected species-level microbial OTUs may represent novel species.Our finding that the proportion of unassigned sequences in proglacial samples was greater than that of lake samples substantiates the significant knowledge gap that currently exists regarding Tibetan glacial biodiversity.
The Tibetan Plateau has warmed 0.3 C/10 years over the past three decades, which is twice the global average warming rate over this period. 63,64Climate warming has considerably increased terrestrial primary production, accelerated deglaciation, and changed regional hydrological regimes on the Tibetan Plateau, 65,66 all of which inevitably impact both terrestrial and aquatic biological communities and ecosystem processes.Glacier-fed waters and their associated biomes are sentinel systems that indicate early changes in the plateau ecosystem.For instance, diatom assemblage composition serves as a sensitive bioindicator of ecosystem integrity in glacial streams and lakes, 67,68 which is also demonstrated in our IDP analysis (Figure 2).Further, warming can induce structural reforms in food webs through trophic cascades and other biotic interactions and ultimately impair glacial aquatic ecosystem functioning. 69,70Therefore, close monitoring of biological communities in glacial meltwater can be a sensitive tool to trace global change effects and gauge mitigation strategies. 69,71ur results shed light on the spatial succession of biological communities, from simple to complex, along glacial flowpaths and provide information regarding the environmental tolerances and requirements of different taxa.The highly heterogeneous community composition across glacier-fed water systems highlights the complexity and dynamic nature of glacier-influenced ecosystems, which remain largely uncharacterized for most glaciers on the Tibetan Plateau.A more extensive approach to profiling the plateau's biodiversity is crucial for gaining a deeper understanding of its organization and dynamics at the community level.Moreover, further analyses linking multiple taxa to environmental alterations and climate warming across glacier-influenced habitats will provide essential information about the effects of global change on plateau biodiversity and aid prediction of future trends of glacial ecosystem processes, stability, and functions.

Limitations of the study
Glacier-fed water systems are highly dynamic ecosystems that show significant seasonal fluctuations.Our study only investigated a single time point during the melt season, while more comprehensive sampling spanning annual variations in community composition will provide information on the temporal effect of glacial influence on the associated aquatic biota.Furthermore, the low species resolution of a large proportion of the organisms detected in the proglacial habitat hampers a thorough knowledge of the taxonomic identities and biological traits of the biota.More research efforts combining morphotaxonomy and DNA profiling would be necessary to fully understand the biodiversity in these understudied ecosystems.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Bioinformatics processing and taxonomic assignment
Each sequencing library underwent distinct bioinformatics processing using the OBITools3 package (https://metabarcoding.org/obitools3) 91 to process the raw Illumina-sequencing reads following a series of sequence filtering steps.Specifically, alignpairedend was employed to align and assemble the forward and reverse sequencing reads, after which sequences with low overlap alignment scores (%0.6) were eliminated using grep.Demultiplexing was achieved using ngsfilter to assign each sequence to its corresponding sample, while obiuniq was used to dereplicate identical sequences.Only sequences with a total count R10 and a length longer or equal to the minimum length (CYA: 350 bp; 708F: 250 bp; BF: 300 bp; Tele02: 120 bp) were retained using obigrep.Finally, obiclean was applied to detect and eliminate sequences that were potentially generated due to sequencing errors or PCR amplification artifacts (r = 0.5).Putative chimeric sequences were eliminated using the uchime3_denovo command in VSEARCH. 72The remaining sequences were then clustered into OTUs at a 97% similarity threshold through the application of sumaclust (http://metabarcoding.org/sumaclust).
We used the assignTaxonomy function in the R package DADA2 73 to assign the cyanobacteria and diatom taxonomies.The SILVA (release 138) reference database was used to classify cyanobacterial sequences, and only OTUs belonging to the phylum Cyanophyta were retained.Both chloroplast OTUs and cyanobacterial OTUs lacking taxonomic information at the class level were also removed.The Diat.barcode (v.10) database 92 was used for diatom taxonomic assignments.The diatom OTUs were initially filtered to include only those belonging to the phylum Bacillariophyta, and only those OTUs with class-level taxonomic information were subsequently retained.
Sequence alignment for invertebrates was conducted using MIDORI ref.
2 (http://www.reference-midori.info/index.html)against the GenBank eukaryotic mitochondrial COI database (release 251). 93All OTUs belonging to the metazoan kingdom (Kingdom Animalia) were initially retained, and then those classified under the subphylum Vertebrata or lacking phylum information were excluded.To enhance identification accuracy, the remaining OTUs were re-matched to the BOLD database (http://www.boldsystems.org/)using the JAMP method (https://github.com/VascoElbrecht/JAMP)with the Boldigger package, 94 as described previously. 95The subsequent results were combined with the Naive Bayesian classifier results to obtain more detailed taxonomic information.
For vertebrate sequence alignment, we utilized MIDORI ref.
2 against the GenBank eukaryotic mitochondrial 12S (srRNA) database.Only OTUs belonging to the Vertebrata subphylum and possessing taxonomic information at the order level were retained for further analysis, and OTUs identified as human sequences were removed from the dataset.To enhance identification accuracy, the remaining OTUs were rematched to the NCBI nonredundant sequence database using the BLASTn program (https://blast.ncbi.nlm.nih.gov/Blast.cgi).If a query sequence matched a single locally occurring species in the database with R98% identity and 100% coverage, that species was assigned to the corresponding OTU.In cases where a query matched multiple local species with R98% identity and 100% coverage, the lowest taxonomic level encompassing all matched species was assigned to the OTU.In addition, OTUs with the same taxonomic name were collapsed.Local species records were accrued from past survey records for fishes, 96 birds, 97 and mammals, 98 as well as available records from the Global Biodiversity Information Facility (https://www.gbif.org/).The sequences of the domesticated pig and dog could not be differentiated from those of their respective wild ancestors, the wild boar (Sus scrofa) and gray wolf (Canis lupus).As wild boars and wolves occur in the study area, we tentatively assigned the sequences to these wild species rather than their domesticated counterparts.
After the taxonomic assignments, the highest count of each OTU in the negative control PCRs of a given library was subtracted from the eDNA PCR results to remove potential contamination.Furthermore, to mitigate the potential impact of tag jumping, any OTU with a read count <10 in each PCR product was eliminated. 99

QUANTIFICATION AND STATISTICAL ANALYSIS
To assess the impact of sequencing depth on the observed OTU richness within each group, we used the rarecurve function in the R package vegan v.2.6-2 to construct rarefaction curves. 74PCR results that had low sequence counts in each group (CYA: <300; 708F: <5,000; BF: <50; Tele02: <5,000) were excluded from subsequent analysis.The predetermined thresholds for inclusion were based on the rarefaction curves (Figure S1) and were consistent with a previous report. 12After all sequence filtering was complete, 47 PCR results were retained for cyanobacteria, 39 for diatoms, 33 for invertebrates, and 28 for vertebrates.Next, the OTU data were normalized using the total sum normalization method 100 to calculate each PCR's RRA, which is the number of reads for a specific OTU divided by the total number of reads of the PCR, averaged across the PCR replicates of each sample.
Diatom-based biological indices are widely employed to assess water quality changes in freshwater environments.The IDP is an effective diatom-based indicator of organic pollution and eutrophication, and it ranges from 0 (low levels of nutrients and organic enrichment) to 4 (high concentrations of organic matter). 29We calculated the IDP using the R package DiaThor v.0.1.0. 75ompositional variation between sites (i.e., b diversity) for different biological groups was quantified using the Bray-Curtis dissimilarity index, which was estimated using the RRA data with the vegdist function in vegan.To assess community compositional differences, we conducted PERMANOVA using the adonis2 function in vegan with 9,999 permutations.The results were visualized using the vegan metaMDS function, which employs NMDS analysis based on Bray-Curtis dissimilarity.
To estimate the phylogenetic relationships among the OTUs of each of the four biological groups, we first conducted multiple sequence alignment of representative OTU sequences using the Muscle method in MEGA X. 76 Subsequently, IQ-TREE2 77 was used to infer the maximum-likelihood tree, and the best substitution model was selected using the ''Auto'' parameter.The tree was visualized and annotated using the Interactive Tree Of Life webtool. 78To gain insights into the specific OTUs that contributed the most to the observed dissimilarities between the proglacial and lake communities within each biological group, we used the simper function in vegan to conduct a SIMPER analysis.
The b-diversity for different biological groups was partitioned into species turnover and nestedness components using the beta.multifunction for Sørensen (qualitative) dissimilarity and the beta.multi.abundfunction for Bray-Curtis (quantitative) dissimilarity in the R package betapart v.1.5.6. 79o identify factors that may influence biological communities, Mantel tests were employed to assess the correlations between biotic assemblage compositional similarity and environmental factors across different sites.The environmental factors encompassed three spatial variables (latitude, longitude, and elevation) and 10 physiochemical water parameters (Table S1).The spatial and environmental data were logarithmically transformed to mitigate skewness, and the RRA data underwent Hellinger transformation prior to analysis.A heatmap was generated in the R package linkET v.0.0.5 97 to visualize the correlations between assemblage compositions and environmental factors.
Both deterministic (i.e., niche-based, such as environmental filtering and biological interactions) and stochastic (including birth/death, dispersal, speciation, and ecological drift) processes can influence community diversity and dynamics, and determining the relative importance of these processes is a central theme in community ecology. 53,101,102To investigate the community assembly mechanisms of different biological groups, a phylogenetic-bin-based null model was applied with the R package iCAMP v.1.5.12. 81By using iCAMP, the relative importance of five assembly processes (homogeneous selection, heterogeneous selection, dispersal limitation, homogenizing dispersal, and drift and others) was quantified based on all pairwise comparisons between samples. 103Phylogenetic compositional variation among assemblages is expected to be low under homogeneous selection (i.e., selection under similar abiotic and biotic conditions), whereas heterogeneous selection leads to a high level of phylogenetic compositional variation.Similarly, homogenizing dispersal refers to the situation where a high dispersal rate leads to low taxonomic compositional variation among assemblages, while dispersal limitation (low dispersal rate) increases community taxonomic variation.Processes other than strong selection and dispersal are collectively designated as drift. 81Therefore, by analyzing the taxonomic (OTU) and phylogenetic (sequence) compositional variations among samples, the relative contributions of different assembly processes can be assessed.The iCAMP analysis was performed using the RRA data and sequences of the detected OTUs and the icamp.bigfunction with a general threshold of 0.2 for phylogenetic signal distance, a randomization number of 1,000, and a minimal bin size determined by the ps.bin function.We conducted this analysis for proglacial and lake samples separately, as these two habitats each showed distinct community compositions for each biological group (see results).Because only cyanobacteria were detected at all proglacial sites, while the other groups were detected at only 0-3 proglacial sites, we assessed proglacial habitat community assembly only for cyanobacteria.Vertebrates were not considered in this analysis because their community assemblies may have been subjected to other variables unfit for iCAMP assumptions.

Figure 1 .
Figure 1.Maps and photographs of the Parlung No. 4 glacier study area and sampling sites (red circles) on the Tibetan Plateau (A) Proglacial sites (PL01-PL05) at the glacier's terminus and along the glacier-fed stream.(B) Lake sites (RW01-RW11) along the glacier's downstream lake, Ranwu.The geological maps were generated using QGIS v.3.22 (https://qgis.org) with an ESRI World Imagery base layer.The insets for (A) and (B) show the geographic placement of the glacier and lake, and the red star indicates the study area.Arrows denote the direction of flow.PL, Parlung No. 4 glacier; RW, Ranwu Lake.Photos show the environments of representative sampling sites.

Figure 2 .
Figure 2. Operational taxonomic unit (OTU) richness (i.e., number of OTUs) distributions of the four biological groups detected at the proglacial (PL01-PL05) and lake (RW01-RW11) sites Vertebrate OTUs included all detected vertebrates, and the richness of domesticated species (cattle, yak, horse, sheep, and chicken) was plotted separately from wild taxa as an indicator of human activity.The scale is log-transformed to facilitate visualization of lower values.The Pampean Diatom Index (IDP) reflects the level of organic pollution and eutrophication.For IDP, light green, light blue, and brown indicate low, moderate, and high nutrient levels, respectively.

Figure 3 .
Figure3.Relative sequence abundance of identified operational taxonomic units (OTUs) for four biological groups along the proglacial headwaters and lake continuum The taxonomic levels displayed are based on more than 70% of the OTUs identified for each group having taxonomic information at that level.The figures depict the 10 most abundant families of cyanobacteria, diatoms, and vertebrates, and the 10 most abundant orders of invertebrates.Low abundance families/orders were collapsed into ''Others.''The OTUs that could not be assigned to a single family (for cyanobacteria and diatoms) or order (for invertebrates) were grouped and are shown as ''Unassigned.''PCR results for certain proglacial sites were excluded due to low read abundance after sequence filtering, thus explaining the blank sections on the plots.

Figure 4 .
Figure 4. Bray-Curtis dissimilarity-based non-metric multidimensional scaling (NMDS) plots of the four biological groups detected in different habitats Blue circles indicate samples from proglacial sites.Green triangles indicate samples from lake sites.

Figure 5 .
Figure 5. Phylogenetic trees in the pie chart are based on 642 OTUs detected for four biological groupsThe inner bar graphs represent the OTU relative read abundance (RRA) in the proglacial (blue) and lake (green) samples.The outer ring indicates the mean OTU RRA across all sites, and the names of the OTUs that contributed R5% to differences between the glacier and lake assemblages are shown in the color corresponding to their more prevalent habitat.

Figure 6 .
Figure 6.Heatmap of correlations (Mantel tests) between the four biological groups' assemblage compositions and individual environmental factors The width and color of the edges represent the Mantel's r value and statistical significance, respectively.Pairwise Pearson correlations between environmental factors are indicated with a color gradient.Lon, longitude; Lat, latitude; Ele, elevation; DO, dissolved oxygen; Tem, water temperature; DOC, dissolved organic carbon; DON, dissolved organic nitrogen; NH 4 + -N, ammonium nitrogen; NO 3 À -N, nitrate nitrogen; TC, total carbon; TN, total nitrogen; TP, total phosphorus.

Figure 7 .
Figure 7. Relative importance of different community assembly processes for cyanobacteria, diatoms, and invertebrates in proglacial and lake habitats The analysis focused solely on cyanobacteria in the proglacial habitat because the sample availability for the other groups was limited.

TABLE
d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d METHOD DETAILS B Study area and sampling B Environmental DNA laboratory procedures B Bioinformatics processing and taxonomic assignment d QUANTIFICATION AND STATISTICAL ANALYSIS