The samples were obtained from the Uraim River in September 2022, which flows through Paragominas. Five sampling points were selected based on varying levels of anthropogenic impact: P1 (P1A, P1B, P1C), P2 (P2A, P2B, P2C), P3 (P3A, P3B, P3C), P4 (P4A, P4B, P4C), and P5 (P5A, P5B, P5C), were collected 1 m away from the banks and in the center of the river for each point. Three biological replicates were collected at each point, and fifteen samples were obtained from all points.

The point P1 was located in an untouched, pristine site, while the final point (P5) was located downstream after passing through the city. The remaining points were intermediate between these two extremes. The samples were collected at a depth of 20 cm, with a volume of 5L taken at each point in triplicate, with help from a Van Dorn bottle. They were preserved on ice and immediately transported to the laboratory. Upon arrival, the water samples were filtered using a negative pressure filtration system with the assistance of a vacuum pump and 47 mm membranes with a porosity of 0.22 µm. These membranes were preserved in 50 mL falcon tubes containing TSE buffer (Tris HCl 10 mM pH = 8, Sodium Chloride 100 mM, EDTA 1 mM pH = 8) and stored at −20°C.

Physical-chemical component measurements were conducted in situ using a multi-parameter water analyzer (HI 9829, HANNA), considering the parameters of pH, electrical conductivity (EC), dissolved oxygen (DO), and total dissolved solids (TDS).

The methodology for mapping the anthropic signature of land use and land cover in rivers employs the concept of a hydrographic basin and geoprocessing techniques. The deforestation level and population density were evaluated following the methodology of Cavalcante et al. (2023). The deforested area was calculated as the sum of all areas with non-natural land use and land cover based on the classification of the Mapbiomas Brasil project (https://brasil.mapbiomas.org/en/) version six for the year 2020.

The sampling points were chosen following a bioprospecting approach, with consideration given to the gradient of anthropogenic impact and accessibility for sample collection. These points were situated at the outset of the pristine site, specifically at coordinates 3° 3′28.512″S, 47° 25′2.172″W. Following the initial point, the subsequent sampling locations were designated as P2 (3° 0′24.737″S, 47° 22′52.941″W), P3 (2° 59′56.264″S, 47° 22′17.735″W), P4 (2° 58′56.895″S, 47° 21′37.938″W) and extending to the river outlet from the city of Paragominas at P5 (2° 57′54.169″S, 47° 20′15.140″W).

In the laboratory, the membranes containing the samples were first shaken at 250 RPM at room temperature. Subsequently, the liquid containing the suspended particles retained in the membrane was centrifuged at 13200 RPM for 10 min in Eppendorf tubes. Genetic material was then extracted following the manufacturer’s instructions of the Pro-DNeasy PowerSoil kit (Qiagen). The integrity of the obtained genetic material was assessed through horizontal electrophoresis. DNA quantification was conducted using Qubit (Thermo Fisher Scientific), and purity indices were obtained via Nanodrop (Thermo Fisher Scientific).

For the construction of DNA libraries, the samples were processed according to the manufacturer’s protocol for the Nextera XT DNA Library kit (Illumina). Sequencing was performed using the NextSeq 550 System High-Output kit, generating paired-end reads of 150 bp, using the Illumina NovaSeq 500/550 High-Output platform.

The quality of the sequencing reads was initially analyzed using the software FASTQC–Version 0.11.9 (Wingett and Andrew, 2018). Subsequently, the reads were trimmed and filtered with a minimum Phred quality score of 20 using the software Trimmomatic (Bolger et al., 2014). The obtained sequences were submitted to the National Center for Biotechnology Information NCBI database (https://www.ncbi.nlm.nih.gov/), with Submission ID SUB14513086 and BioProject ID PRJNA1122411, accessible with the following link https://www.ncbi.nlm.nih.gov/sra/PRJNA1122411.

The pathogenic species in the collections were identified by comparing the output from Kraken2 with a local database of pathogenic bacteria from various isolates associated with pathogenicity in food, water, and clinical hosts. Analyses were conducted at both the species and genus levels. The local database was created using complete genome records labelled as pathogenic by the NCBI Pathogen Detection Project (https://www.ncbi.nlm.nih.gov/pathogens/). Records were extracted using the NCBI FTP protocol and manipulated using the Python programming language. For information extraction from HTML pages and for the grouping, filtering, and sorting of records, the beautifulSoup and pandas packages were utilized, respectively. The records were saved in the MongoDB database, employing the NoSQL format, hosted in a Docker container https://github.com/labgm/pathogen-detection.git.

For the identification of different taxonomic groups, the previously trimmed reads were used. These reads were then evaluated using the Plants, Fungi, and Protists (PFP) database within the Kraken2 program (Wood et al., 2019). The output format of Kraken2 was converted to be recognized by the MicrobiomeAnalyst program (Chong et al., 2020) for the transformation of abundances using Trimmed Mean of M-values (TMM), then, relative abundances were processed to determine alpha and beta diversity indices and associations between taxa.

For the identification of bacterial virulence factors, the Virulence Factor Database (VFDB) (Liu et al., 2022) was used. This database contains 32 genera of bacterial pathogens, classified into 14 basal categories. The database was indexed using the Bowtie2 software (Langmead and Salzberg, 2012), and then the reads obtained from the sequencing were aligned against the database to capture the presence and abundance of virulence factors in each sample. The alignment result was processed using the Pileup software to identify and quantify the aligned reads at each position of the reference database. The information obtained from Pileup, including the columns of identification and abundance, served as the basis for the subsequent step where they were organized into a tabular structure file containing the identified virulence factors and their respective abundance values for each sample. This facilitated the comparison and interpretation of the data.

The data obtained from the different analyses were evaluated for normality using the Shapiro-Wilk test. Subsequently, they were analyzed using ANOVA or Kruskal–Wallis tests, multiple regression, and principal component analysis (PCA) using the PAST 4.05 software (Hammer et al., 2001). For the preparation of the graphs, the Orange Data Mining program (Demsar et al., 2013), as well as the Seaborn and Vegan-Altair libraries of Python, were used.

The points where samples were collected in the Uraim River are highlighted in yellow in Figure 1, observed through satellite images. In Figure 1A, the population density in inhabitants per square kilometer along the river is emphasized, allowing for the observation of human settlement and agricultural areas near certain areas of the river. In Figure 1B, the percentage of deforestation in the areas along the river is further observed. It is noted that the collection sites showed a deforestation percentage ranging from 70% to 73%, with population densities ranging from 1 to 184 inhabitants per square kilometer.

The Uraim River basin exhibits considerable deforestation of its native vegetation, resulting in a 25% reduction in average precipitation during the deforested period compared to years with forest presence, consequently leading to a 22% reduction in channel flow (Sardinha and Ventura, 2022). Different anthropogenic activities have repercussions on various water resources (Pereira Júnior et al., 2023a).

In the municipality of Paragominas-PA, the forest area was 3420.27 km² in 2000 and decreased to 2784.23 km² in 2021, representing a forest area loss of 18.6%. Regarding the expansion of the urban area, the development began in 2000 (133156 km²) and 2021 (20.2687 km²), which is linked to population growth: in 2000, the population of Paragominas was 76450 inhabitants, and it increased to 115538 inhabitants in 2021. This implies changes in ecosystems, such as the reduction of the precipitation rate that affects the water recharge of rivers, streams, and creeks, among others (Pereira Júnior et al., 2023a). These characteristics are evident in the images obtained for this study and also in the alterations in chemical and microbiological components.

In the physicochemical analyses, acidic pH was found, ranging between 3.547 and 4.627, with significant differences between the points (ANOVA; F = 23.47; p-value = 0.0001). notably, these differences occurred in the comparison between: points P1 and P2 (Tukey’s test; Q = 9.235; p-value = 0.001); P2 and P3 (Tukey’s test; Q = 6.527; p-value = 0.011); P1 and P4 (Tukey’s test; Q = 7.707; p-value = 0.004); P3 and P4 (Tukey’s test; Q = 4.999; p-value = 0.048); P2 and P5 (Tukey’s test; Q = 2.014; p-value = 0.0003), and P4 between P5 (Tukey’s test; Q = 9.721; p-value = 0.001) (Supplementary Figure S1).

The values of electrical conductivity varied across all points but were not significant. In the case of dissolved oxygen, a decrease was evidenced at P5, with average values of 4.413, showing significant differences (ANOVA; F = 8.148; p-value = 0.006) for this point compared to the others. For total dissolved solids (TDS), despite the variation in values at the collection points, the differences were not significant (Table 1).

These findings underscore a potential public health concern should the population consume untreated water, including water from the river, typically used for domestic, agricultural, and livestock purposes. Anthropogenic activities stand out as primary factors influencing the quality of surface waters (Tornevi et al., 2014; Vinueza et al., 2021). In the Uram River area, agricultural practices have been observed to increase the concentrations of various chemical compounds in the river. As it flows through urban areas, the river’s natural pH, nitrogen and oxygen concentrations, turbidity, and other parameters can be altered, affecting the standards of freshwater quality intended for human consumption (Pereira Júnior et al., 2023b).

In the sequencing, an average of 35468174.13 reads were obtained from the different sample collection points, of which about 33% were used for taxonomic classification. Bacteria accounted for 30% of the reads, while opportunistic pathogenic bacteria for humans represented 0.34%, Table 1.

Regarding total diversity in these ecosystems, up to 7.190 different species were identified, with Shannon index values above 7.184 and Simpson index values close to 1, indicating no dominance of a single species (Table 1). In the Tama River in Tokyo, Japan, it was observed that bacterial communities exhibited significant differences along the river’s course, showcasing variations in diversity (Mizusawa et al., 2021). This phenomenon was also noted in the Uraim River in Paragominas, Brazil, highlighting the dynamic nature of microbial communities in urban river systems. Interestingly, the Amazon region, home to the Amazon River, lacks comprehensive studies on microbial diversity. Despite being the world’s largest river in terms of discharge, the Amazon River remains relatively understudied, with existing research exhibiting considerable heterogeneity (Li et al., 2021).

In our study of the Uraim River, the most abundant phyla were Pseudomonadota and Actinomycetota, the taxonomic class of Betaproteobacteria was noteworthy, followed by Alphaproteobacteria, Actinomycetes and Gammaproteobacteria. However, in freshwater studies, the dominant classes are usually Alphaproteobacteria, Betaproteobacteria, and Actinobacteria, which complex changes in nutrient composition can alter (Santos-Júnior et al., 2020; Mizusawa et al., 2021).

A higher abundance of pathogenic species was detected at point P5, averaging 226.7 (Figure 2A). The Shannon index showed lower values on average at point P2, with a value of 3.955 ± 0.592, and the highest average value of 4.312 ± 0.102 was found at point P3 (Figure 2B). The Simpson index presented average values between 0.914 and 0.971 at the same aforementioned points (Figure 2C).

Beta diversity analysis of biological replicates revealed a clustering trend in the bacterial composition corresponding to P5, with R ² values of 0.576 (PERMANOVA F-value = 3.3961; p-value = 0.001), indicating a cumulative variation of 70.6% for the PCoA analysis (Supplementary Figure S2).

Regarding the phylum composition of pathogenic microorganisms identified across all sampling points, Pseudomonadota, Bacillota, and Campylobacterota are highlighted, respectively. Among the eight microorganisms identified, Gammaproteobacteria were predominant, followed by Betaproteobacteria and Bacilli, with the remaining classes taxonomic in smaller proportions (refer to Supplementary Figure S3). Among the 27 identified families, Enterobacteriaceae, Burkholderiaceae, and Pseudomonadaceae exhibited the highest abundance (see Figure 3A). Particularly for Enterobacteriaceae, a substantial increase was observed in points 4 and 5, with statistically significant differences noted between these points (ANOVA; F = 78.96; p-value = 1.59 × 10⁻⁷) (see Figure 3B).

Metagenomic analysis allows for the simultaneous screening of various potential pathogens and fecal indicators. The river’s microbiome serves as a significant indicator of urban pollution alteration, demonstrating a prevalence of human intestinal bacteria, predominantly pathogenic (Acharya et al., 2022). Among the most abundant genera identified at all detecting points were Burkholderia and Pseudomonas, with a notable presence of Acinetobacter, Escherichia, Klebsiella, and Enterobacter at point 5 (Figure 4A). The most abundant pathogenic species belong to the ESKAPE group. Less frequently detected species include Burkholderia cepacia complex (Bcc), E. coli, and Stenotrophomonas maltophila (Figure 4B). Regarding ESKAPE group species, Enterobacter spp. Exhibited a distribution with minor variations across sampling points (R ² = 0.768; p-value = 1.825 × 10⁻⁵), consistent with the observations for Acinetobacter baumanii and Klebsiella pneumoniae, which displayed higher abundance at point 5. Notably, Pseudomonas aeruginosa demonstrated an inverse trend compared to the species above (Figure 4C), which could influence changes environmental characteristics and microbial communities (Philippot et al., 2021), but we could not identify specific facts against these microorganisms due to the complexity of the variables that make up the ecosystems.

In the correlation analysis with the relative abundances of pathogenic bacterial species, Enterobacter cloacae stands out as significantly associated with the presence of Klebsiella variicola, Citrobacter freundii, Aeromonas veronii and Enterobacter kobei. The latter is correlated with the presence of E. coli and is involved with the occurrence of Acinetobacter nosocomiais, Streptococcus susis, Enterococcus hirae and Enterococcus faecium, which are associated with Streptococcus pyogenes and Streptococcus equi, as indicated by the interaction network presented in Figure 5.

The substantial presence of human pathogens is notable, constituting approximately 7% of all taxonomically described bacterial species. These pathogenic entities span 24 distinct taxonomic classes, with prominent contributions from Gammaproteobacteria and Actinomycetia in terms of genus and species diversity. The identified pathogen pool encompasses 1513 species categorized into 327 genera, of which 267 genera encompass at least one established pathogenic species (Bartlett et al., 2022). Waterborne illnesses persist as a significant global health concern, prompting numerous regulatory bodies to utilize coliform assessments in water quality testing as an indicator of fecal contamination levels and associated disease risks (Hamner et al., 2019).

The findings presented in this study align with the expected bacterial populations that thrive in sewage-contaminated environments, characterized by a dominance of Proteobacteria and Bacteroidetes (Godoy et al., 2020). Additionally, the prevalence of Actinobacteria in freshwater environments is consistent with previous studies (Mizusawa et al., 2021).

The presence of pathogenic microorganisms allows us to assess the quality of these aquatic systems in terms of environmental impact and anthropogenic activities. In this study, the presence of the ESKAPE group is notable, given its ability to possess and transfer antibiotic resistance genes (Denissen et al., 2022; Aguilar-Salazar et al., 2023). The results highlight a higher abundance, on the bodies of water, of K. pneumoniae, E. cloacae, E. faecium, and E. faecalis. Additionally, the metagenomic analysis reveals the presence of bacteria that serve as bioindicators for chronic diseases in humans, such as Fusobacterium nucleatum associated with colorectal cancer (VanEvery et al., 2023).

Upon examination of the virulence genes identified across the various sampling points, as illustrated in Figure 6, it was noted that point P5 exhibited a greater prevalence of unique genes, not identified in other points. This observation aligns with expectations, given that P5 is situated downstream from the urbanized area of Paragominas, experiencing a range of anthropogenic influences and cumulative pressures from upstream points. Additionally, a total of 315 genes were found to be common among all five sampling locations, indicating a potential presence of native core microbiota unaffected by anthropogenic activities.

Virulence factors (VFs) can affect a wide range of cellular processes, such as cell-cell signaling, ion secretion, protein synthesis, mitosis, cytoskeleton structure, and mitochondrial function, although their role in natural environments is not yet well understood (Vinueza et al., 2021). Our analyses identified numerous genes related to virulence factors, especially in sites with higher anthropogenic activity (Figure 6), particularly in P5 where there is a greater abundance and diversity of VFs. This site also showed a higher quantity of pathogenic bacteria (Supplementary Figures S4A, B).

The diversity of virulence genes calculated was highest at point P5 with a value of 105.6 for the Fisher alpha diversity index, followed by 77.96 at point P4, values close to 60 at points P1 and P2, and the lowest value of 58 at P3 (Figure 7A). The most relevant virulence mechanism is adhesion, present in all points; points P2 and P5 highlight the effector delivery system mechanism, and in P1, the nutritional/metabolic factor and invasion mechanisms are present (Figure 7B).

In studies conducted on rivers traversing urban areas, it has been revealed that urbanization increases environmental risk and the occurrence of resistance and virulence genes compared to those in rural regions. Among the main virulence factors (VFs) are offensive, nonspecific, defensive, and regulatory VFs, which were found in all samples. The dominant subtypes include iron acquisition systems, adhesion, secretion systems, regulation mechanisms, toxins, and antiphagocytosis (Mizusawa et al., 2021; Li et al., 2022). In our samples, secretion systems and adhesion were the most relevant factors.

In all sampling points, a high abundance of the gene fbpC, which is involved in the Secreted fibronectin-binding protein C (VF0311), was observed. These adherence proteins can ensure the adhesion of these pathogens to hosts, with many of these molecules being diverse and not fully described (Aono et al., 2023).

At sampling point P5, there is a notable abundance of genes involved in secretion systems, with the most prominent being the gene vgrG/tssI, which is associated with the synthesis of type VI secretion system tip protein VgrG (VF1337). At P2, the gene vgrG1a, which synthesizes the type VI secretion system substrate VgrG1 (VF0334), stands out (Supplementary Figure S4A).

The identification of genes related to Type III and Type IV secretion systems is essential in bacterial pathogenesis, facilitating direct interactions with host cells. The presence of these genes may reflect regional variations in the bacteria-host dynamics or environmental factors that favor certain pathogenic strategies (Rout et al., 2023). Metagenomic analysis of the Ganges River (India) revealed the widespread distribution of genes associated with flagella, which are related to Type III and Type IV secretion systems, the latter being the most abundant marker identified in different sites. This study highlights the complexity of bacterial virulence and its dependence on environmental context. The detection of flagella-associated genes in various locations implies a widespread reliance on motility and adhesion as critical virulence factors. This suggests a common strategy employed by bacteria to establish infection and colonization across different geographic areas (Godoy et al., 2020; Rout et al., 2023).

In point 1, the gene lic2A synthesizes glycosyltransferase (VF0044), and glnA1 is involved in the production of Probable glutamine synthetase, type I (VF0816) (Supplementary Figure S4A). The glnA1 gene is involved in nitrogen metabolism and is necessary for the development of Mycobacterium tuberculosis and E. coli in the environment (Xu et al., 2023).

The secretion systems identified in this study are characteristic of a variety of Gram-negative bacteria, contributing to the virulence of Type I and Type VI secretion systems, which are commonly described in the genome of P. aeruginosa (Green and Mecsas, 2016), and in members of the genus Pseudomonas, Acinetobacter, Klebsiella, and Enterobacter, which are part of the ESKAPE group and they are the most abundant among the identified pathogenic bacteria identified in our results, which could contribute to the high abundance of this gene.

Another mechanism of adaptation and resistance is the production of biofilms, whose production is induced by the presence of N-acyl-homoserine lactone (Alio et al., 2023), providing significant advantages to bacteria such as resistance to antibiotics and protection against host immune responses. In the environment, biofilms allow bacteria to resist unfavorable conditions (Gomes et al., 2020; Rout et al., 2023), and are present across different sites studied (P1, P2, P3, P4 e P5) in Paragominas.

In the PCA analysis, an accumulated percentage of 88.04% of variation in the variables was found, with 58.8% variation in principal component 1 (PC1), highlighting the following associations: between the physicochemical components of pH, Ec, and TDS with the alpha diversity indices of Shannon, although these were not evident in one analyzed point. In the case of P5, a strong correlation was observed with population occupancy, and Fisher’s alpha diversity of virulence genes (Pearson; ρ = 0.972; p-value = 0.006). Furthermore, a high correlation was identified between the relative abundance of the ESKAPE group and the percentage of deforestation (Pearson; ρ = 0.86; p-value = 0.061); a significant inverse association between dissolved oxygen variables and deforestation percentage (Pearson; ρ = −0.912; p-value = 0.03); and a significant association between virulence gene diversity and the percentage of inhabitants per surface area (Pearson; ρ = 0.972; p-value = 0.005) and with deforestation (Pearson; ρ = 0.907; p-value = 0.033) (Figure 8).

The significant associations found between urban and rural population growth and the functioning and biodiversity of the ecosystem by microorganisms responsible for the purification and removal of nutrients in rivers, as well as the nutrient cycle and interactions between plants in the soil (Gomes et al., 2020), are reflected in different analytes. Similarly, biochemical and chemical oxygen demands are influenced by the urbanization process, temperature, dissolved oxygen, ammonia, and phosphorus (Giuliatti et al., 2017).

In this study, the assessed physicochemical properties underwent significant alterations, acting as crucial factors that modify bacterial community structures within these ecosystems and facilitate the acquisition of genetic components associated with their adaptation, including virulence genes for higher organisms (Godoy et al., 2020). Among the identified components, DO is inversely associated with the diversity of virulence genes and the presence of certain pathogens. This association could be attributed to anthropogenic activities, with human presence notably increasing the likelihood of greater contact with these types of organisms.