Prokaryote communities along a source-to-estuary river continuum in the Brazilian Atlantic Forest

Study area

Three sites were chosen along the Juliana River in the southeastern part of Bahia State, Brazil. The river drains the most important watershed in the region in terms of size and economic and ecological significance. Currently, the Juliana River is located entirely within a legally protected area, the Environmental Protection Area of Pratigi (Fig. 1). Its basin comprises an area of 299.8 km², through which the river runs almost linearly over 47 km. The source is in the Papuã Mountains. Several tributaries join the river along its way to the Serinhaém estuary (Mascarenhas et al., 2019; Ditt et al., 2013), where the city of Ituberá is located, a small urban area with less than 30,000 people where tourism is the main economic activity (IBGE–instituto brasileiro de geografia e estatística, 2020). Ituberá has been constructed within a mangrove forest, which has been retained along urban waterways and mudflats (de Santana et al., 2021b). In contrast, most of the upstream reaches enable the observation of minimally impacted environments, because the upper portions of the watershed are considered to be highly conserved, lending themselves to ecological, hydrological and biogeochemical research. This includes studies of the biodiversity and ecology of microbial communities in mangrove sediments (de Santana et al., 2021a).

The Juliana River basin is subdivided into three administrative sections, I, II and III. Section I corresponds to the highlands of the Papuã Mountains. A site located there has been designated the Source site for the purpose of the present study. Section II corresponds to the downstream Valley region, which is mostly dominated by forest cover interspersed with a few agroforestry systems. Section III is the lowermost part of the hydrographic basin, hosting ecosystems ranging from tropical forest fragments to mangroves (Mascarenhas et al., 2019), including in and near Ituberá City close to where the sediments were collected. Nevertheless, this area still experiences little direct impacts by industrial development and family farming predominates land use (da Silva Pereira et al., 2022).

The field study presented here was approved by the state government organization Fundação de Amparo à Pesquisa do Estado da Bahia (project number: FAPESB/CNPq n° 794014/2013; permit number: 794014/2013). Portions of this text were previously published as part of a doctoral thesis (de Santana, 2020).

Sampling and genomic analyses

Sediments were collected in February 2019 at the three sites selected in the Juliana River (Source, Valley, and Mangrove). At each site, three collection points were selected. Each site had to be at least 1.5 m apart from one another, free of visual vegetation, contamination, or pollution, and at a margin of the river where water depth exceeded 10 cm. Sediments were collected with a cylindrical core sampler, taking precautions to avoid disrupting rhizospheres associated with vegetation. Each sample consisted of a 10 cm deep surface sediment sample. Plant litter and other coarse particulate organic matter was manually removed from the core before placing the sediment samples in plastic bags on ice in thermal boxes and immediately transporting them to the laboratory for chemical and genomic analyses.

Physical-chemical parameters such as temperature, pH, conductivity, and dissolved oxygen in the water column were measured at each site using a multiparameter probe (YSI model 85; Yellow Spring Instruments Inc., Yellow Springs, OH, USA). Additional environmental variables such as concentrations of Pb, Zn, Cu and Cd at each site have been previously reported (da Silva Pereira et al., 2022; Mascarenhas et al., 2019; Table S1). Since Cd concentrations were below detection limit at all sites, this variable was not included in the data analysis. In the laboratory, an aliquot of each sediment core was frozen at −20 °C for subsequent DNA extraction, while the remainder of the sample was used to measure organic matter (O.M.) content.

The total genomic DNA was extracted from 0.25 g of sediment using the PowerSoil DNA Isolation Kit (Qiagen, Germantown, MD, USA) and stored at −80 °C before analysis. After DNA extraction, the samples were sent on dry ice to Novogene Bioinformatics Technology Co. Ltd., Beijing, China for amplification of bacterial 16S rRNA genes, using the 515F and 806R primers (Table S2), followed by Illumina NovaSeq 6000 paired-end (2 × 250) sequencing (Thompson et al., 2017). Since sequencing of one of the samples from the Valley site failed, analyses were limited to the two remaining replicates.

Trimmomatic (Bolger, Lohse & Usadel, 2014) was used to filter and trim the demultiplexed sequences (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:100). All reads were subsequently denoised using DADA2 (Callahan et al., 2016) in QIIME2 (Bolyen et al., 2019), merged using QIIME2 (File S1; Table S3), and then clustered into amplicon sequence variants (ASVs) (Table S4). Alpha-rarefaction was calculated using QIIME2 (Fig. S1) and set to 41,000 reads for the purpose of alpha- and beta-diversity analyses (Figs. S2, S3). All diversity analyses were performed using QIIME2’s default parameters (File S1).

Statistical analyses

Taxonomic assignment was performed using QIIME2’s naive Bayes scikit-learn classifier (Bokulich et al., 2018) trained with the 16S rRNA gene sequences in the SILVA database (SILVA 138-99-515-806) (McDonald et al., 2012). The taxonomic feature table (Table S5) was resolved to the genus level for analysis (Table S6) using QIIME2. For each site, a bar chart was made of the phylum and class using the mean percentage of taxa abundance was calculated across replicates (Fig. 2A; File S2). Classes of high relative abundances (2% of the total community per site) and phyla were identified, and a heatmap of relative genus abundances generated for each replicate sample (Fig. S4; File S3).

Taxa resolved to the genus level were considered common across sites if they accounted for at least 0.1% of the reads per site, occurred in at least two replicates per site, or represented at least 1% of the reads in a single replicate. These criteria had to be met for each of the three sites (Fig. 2B; File S4; Table S7).

To determine how many taxa, resolved to the genus level, were only found at any given site, we first required each taxon to be minimally present at only one site. Minimal presence was defined as being greater than 0.001% of the total population per site, or being, on average, greater than 0.0001% of the population per site per replicate (Fig. S5; File S5; Table S8).

A site-specific analysis of significant differential abundances was performed using the ANCOM-BC package in QIIME2 (Tables S9–S11). We further subset these taxa to identify those that were significantly distinct to a single site, relative to the other two sites (ANCOM-BC, q-value < 0.01) and that also represented a substantial percentage of the total population at that site (>1% total population), (File S6; Fig. 3D). A Venn diagram showing the overlap of significantly different taxa at each site is also available as a supplemental figure (Fig. S6).

The Vegan package (Dixon, 2003) was used to test correlations between community structure and environmental variables in R environment (version 4.2.2). Distances were calculated using metaMDS (distance used was Bray-Curtis) (Fig. S7; Table S12; File S7) and environmental variables were fit using envfit (Fig. 4B; Table S13; File S8).

The sequencing data is available from NCBI BioProject PRJNA650560. The entire computational workflow is available in a GitHub repository: https://github.com/pspealman/Project_Juliana_River_basin.

Taxonomic composition of sites and predominant groups

After quality filtering and taxonomic assignment, the 879,453 sequences remaining displayed the following pattern: 91.0% of the reads were associated with the kingdom Bacteria, 8.3% were associated with the Archaea and 0.6% were not assigned to either of these prokaryotic kingdoms. In total, ASVs were assigned to 85 phyla, 202 classes, 457 orders, 699 families, 1,089 genera and 458 species (Table S4).

We identified 18 highly abundant classes with a mean abundance per site of at least 2% (Fig. 2A). These classes constituted nine bacterial and two archaeal phyla. The two archaeal phyla, Crenarchaeota and Thermoplasmatota (as well Halobacterota, which was just below the 2% cutoff) were present at all sites, although they were most frequent in the mangrove sediments. For the Bacteria domain, the three sites shared similar dominant phyla, with Proteobacteria exceeding 10% and Bacteroidota, Bacillota (Firmicutes), Chloroflexota, and Desulfobacterota accounting each for >5% at all sites. Combined, these five phyla and their 11 classes represented the majority of the prokaryotic populations (50–64%) at each site.

This large overlap prompted us to assess how many of the more abundant genera were present at all sites (see Methods). We found 87 such taxa, 77 of which were resolved to the genus level (Table S7; Fig. 2B), which together accounted for 72% (Source) and 61% (Valley and Mangrove), respectively, of the total abundance and could thus constitute the core microbiome in sediments of the river.

Community differences among sites

ANCOM-BC analysis indicated that abundances of numerous taxa significantly differed between pairs of sites (Figs. 3A–3C; File S6). The greatest difference occurred between the Source and Mangrove sites (Fig. S6; Tables S9–S11). Genera specific to only one of the study sites (Fig. S5) included 87 taxa that were unique to the Source site, two to the Valley site, and 63 to the Mangrove site. However, these taxa represent very small proportions of the total communities, with 0.65% being unique to the Source site, 0.03% to the Valley and 1.1% to the Mangrove site (Table S8). Resolved to the genus level, some taxa were significantly more abundant at one site compared to the two others (ANCOM-BC, q-value < 0.01) and represented a notable percentage of the total abundance at that site (>1% total population). We found nine such taxa at the Source site and eight at the Mangrove site (Fig. 3D), whereas none were more abundant at the Valley site, although sediments at that site had more reads that could not be assigned to any taxon (‘Unassigned’).

Community structure, diversity and environmental variables

Prokaryotic diversity expressed as the Shannon entropy index was highest at the Mangrove and lowest at the Valley site (Fig. 4A); however, site differences were only significant in the omnibus test (p = 0.04). Similarly, differences in community composition between sites assessed by the Weighted UniFrac distance measure (Fig. S3) were only significant in the omnibus PERMANOVA (p = 0.007). Site differences among the prokaryotic communities are also shown in the PCA, which separated the Source site from the Valley and Mangrove sites along PC1 (Fig. 4B), with copper (Cu) concentration as the most influential environmental variable (p = 0.011). Nearly significant differences in the concentration of zinc (Zn) (p = 0.063) were primarily related to PC2, whereas temperature, dissolved oxygen, organic matter (O.M.), Ni, salinity, Cr, pH, and Pb had no significant influence.