Rock Surface Colonization by Groundwater Microorganisms in an Aquifer System in Quebec, Canada

Abstract

Aquifers are rich in microbial diversity. However, there is a lack of information about sessile communities in these environments because of the difficulty in sampling fresh in situ rock surfaces. Thus, this study’s objective was to better understand the sessile community in a fractured aquifer. Additionally, the impact of the rock mineral composition on microbial community composition during colonization was explored. Using a system of bioreactors, we recreated the environmental conditions of a 1.5 m deep aquifer in Covey Hill (QC, Canada) using groundwater samples collected from the site. We carried out 16S/18S rRNA amplicon sequencing of the water and sessile communities after 24 days of incubation. Our data showed that many microbial taxa overlapped between the sessile and planktonic communities, indicating colonization of the solid surfaces. Quartz and feldspar had a significant impact on bacterial community structure. Sessile communities were dominated by Gaillonella, Alkanindiges, unclassified Acetobacteraceae, Apoikiales, Glissomonadida, and Synurales. We could not detect any Archaea in the sessile community. The sessile communities contained bacterial genera involved in iron cycling and adapted to acidic and low-carbon-concentration environments. Eukaryotic predators dominated the sessile community.

1. Introduction

The subsurface contains up to 40% of the world’s prokaryotic biomass [1]. These microorganisms play an important role in mediating biogeochemical cycles, such as the carbon cycle and the nitrogen cycle, by recycling inorganic and organic compounds [2,3,4]. Aquifers are regions in the subsurface where groundwater circulates between rocks and solid particles and are an important reservoir of freshwater used by many people around the world [5]. Groundwater flowing in aquifers can be used for human consumption, irrigation, agriculture, and much more [6], which is why it is important to understand how to protect these environments. These dark and nutrient-poor habitats contain microorganisms with two distinct lifestyles: planktonic and sessile [6]. The planktonic microorganisms are free-living in the groundwater, while the sessile microorganisms are attached to solid surfaces. Sessile microbes have the advantage of accessing minerals on rock surfaces, which some species can exploit [7], as well as accessing organic matter that tends to accumulate on surfaces [8]. Nutrients can be released from minerals and can be directly used by microorganisms, while others form precipitates that must solubilize to become available [9]. The dissolved nutrients may have positive or negative effects on the microbial community. Also, as minerals are heterogeneously distributed on rock surfaces, habitat selection occurs, influencing the structure, diversity, and phylogenetic variability of microbial communities in aquifers [9,10].

Different rock types show varying levels of colonization. For example, limestone and dolomite seem to be better colonized by microorganisms than feldspar or shale [11]. Furthermore, microorganisms tend to prefer rocks containing important nutrients like phosphorus and iron [7]. Positively charged rock surfaces can attract the negatively charged bacterial cell wall [7,9,12]. Moreover, metabolic processes carried out by the microorganisms can alter the rocks’ surfaces and thus release elements that can attract chemotactic organisms to the rock [7,13].

Unfortunately, the difficulty in obtaining in situ rock surface samples from aquifers makes it complicated to study these sessile communities, and thus most studies have focused on the planktonic community [4,14,15]. However, previous studies have showed that sessile communities are more diverse than their planktonic counterparts [15,16,17] and are significantly more abundant [18]. Indeed, planktonic bacteria make up only 0.01–10% of the total bacterial community or biomass in aquifers, with the rest being sessile bacteria attached to particles or rock surfaces [19]. Therefore, to better understand assembly and composition of both communities present in a fractured 1.5 m deep aquifer system, we used a set of bioreactors and 16S/18S rRNA amplicon sequencing to compare the planktonic and sessile communities. We also studied the influence of the rock’s mineralogical composition on the microbial colonization of rock surfaces. We hypothesized that the sessile community would comprise microorganisms that have transitioned from a planktonic state, and that the geological characteristics of rocks would impact the microbial community composition.

2. Materials and Methods

For this study, we used a system of three bioreactors (used as replicates) to recreate the environmental conditions and settings of the studied aquifer ecosystem, in the laboratory. Groundwater was collected in the field and was pumped through the systems containing sterile rock chips, to establish colonization of the rock surfaces by groundwater microorganisms. These bioreactors allowed us to circumvent the difficulties of repetitive and invasive sampling of in situ sessile populations in numerous subsurface samples. We monitored dissolved oxygen and pH daily in our in situ experimental setup to ensure the stability of the abiotic conditions. However, it was not possible to completely recreate natural conditions, such as the buffering capacity of the rock layers or continuous slow seeping of water from the surface.

2.1. Study Site, Sampling, and Water Geochemical Measurements

The sampling site is a fractured aquifer (45°00′27.4″ N 73°49′07.5″ W) in Covey Hill, located in the Saint-Lawrence lowlands in the Quebec province of Canada (Figure 1). This superficial aquifer system is composed of deformed and fractured sandstone [20], which outcrops up to the surface, making it an ideal system to study since no drilling is needed to sample aquifer rocks. Furthermore, the aquifer is hydrologically connected to a surface peat bog [21], which is likely the reason why the groundwater is typically acidic [22].

Rock slabs and groundwater were collected in November 2021. The different underground layers of the Covey Hill aquifer system are mostly made of felspar [23], which outcrop up to the surface. Thus, the rock samples at the surface have a similar mineralogical composition compared with those found in the underground aquifer ecosystem. Groundwater samples were collected from a 1.5 m deep well after pumping out a threefold volume to remove stagnant water using 12 Volt Mini-Monsoon Plastic Pump (P-10300, QENEQ, Montreal, QC, Canada). We collected 30 L for the bioreactor incubation experiments (3 times 10 L per bottle) in sterile PYREX low-actinic glass bottles. Also, 1 L was taken four times in sterile polypropylene bottles (Nalgene, Sigma-Aldrich, St. Louis, MO, USA): before collecting the first 10 L to determine the initial planktonic community in the groundwater, and in between the collection of each 10 L sample after that. An extra liter was taken from a nearby peatbog, which is the major recharge source for the aquifer system [21]. All bottles were autoclaved prior to sampling. Groundwater temperature, pH, and percentage of dissolved O₂ (DO) were measured in the field using a multiparameter probe (OAKTON PD 450, Cole-Parmer, Vernon Hills, IL, USA.

Furthermore, the following geochemical parameters were measured in the initial groundwater: dissolved inorganic carbon (DIC), dissolved organic carbon (DOC), NH₄⁺, NO₂⁻, and NO₃⁻. For measurements of DIC and DOC, the water was collected in gas-free glass bottles using 0.2 μm filters. For measurements of the nitrogen compounds, aliquots of water were filtered in the field using either 0.45 μm filters (for the NO₂⁻ and NO₃⁻) or 0.2 µm filters (for NH₄⁺). All the samples were analyzed at the UQAM laboratory of the Interuniversity Research Group in Limnology (GRIL). For DIC and DOC, samples were analyzed with OI Analytical Aurora 1030W TOC Analyzer (Beckman, Fullerton, CA, USA) by using a persulfate oxidation method. For the NO₂⁻ and NO₃⁻ measurements, the samples were analyzed with a continuous flow analyzer (OI Analytical Flow Solution 3100, Beckman, Fullerton, CA, USA) using an alkaline persulfate digestion method coupled with a cadmium reactor. Ammonium measurements were taken with a Flow Solution 3100 autosampler by using a chloramine reaction with salicylate to form indophenol blue dye (EPA Method 350.1). These geochemical parameters were also measured during the incubation experiments in the water, every two days from day 6 of the experiment to day 24.

2.2. Rock Chip Preparation and Mineralogic Composition

The outcropping rock slab collected prior to groundwater sampling was first cleaned with soap and water. Seventy-two rock chips were then prepared by drilling the rock slab into small core plugs (18″ nova voyager dvr, King, Dorval, QC, Canada), and sawing the plugs into rock chips approximately 1.3 cm in diameter and measuring 0.3 cm large (DREMEL 3000, DREMEL, Mount Prospect, IL, USA). The 72 rock chips were further cleaned by vortexing them separately in a HCl 10% solution for 2 min and by sonicating them with a probe (Sonifier Cell Disruptor 185, Branson, Rungis, France) for 1 min. The pellets were then autoclaved, and 24 chips were inserted in each bioreactor (Supplementary Material Figure S1).

To determine the mineralogic composition on each rock surface, pictures of both flat sides of each chip used in the incubation experiments were taken. We used a Leica LED2000 stereo binocular microscope (Leica Microsystems, Wetzlar, Germany) that had a camera attached to it. The pictures were then analyzed using two programs: first Ilastik on the Pixel Classification setting and then ImageJ v.1.54. Identification of the minerals on the rock surfaces was first carried out manually on the basis of physical characteristics such as color, luster, cleavage, and habit, for the program to then be able to automatically classify the rest of the surfaces. Thus, the identification process was faster and more reliable than identifying the minerals ourselves. Using ImageJ, we were able to determine the percentage of each mineral on the surface of the rock using the Threshold tool.

2.3. Bioreactor Incubation Setup

To recreate an environment similar to the aquifer system in the laboratory, we used a bioreactor system (CBR 90 Standard CDC Biofilm Reactor, BioSurface Technologies, Bozeman, MT, USA). Three separate bioreactors were set up, acting as triplicates to check the results’ reproducibility (Figure 2). Each bioreactor contained 8 columns with 3 rock chips; therefore, they contained 24 sterile autoclaved rock chips in total at the start of the experiment. The groundwater collected in the field was pumped at a rate of 0.289 mL/min (IPC, Ismatec, Malente, Germany) from the 10 L field bottles into the bioreactors for 24 days, and the experiment acted as an open system, with the groundwater not being recirculated twice. The excess water was collected in other sterile 10 L containers. Simultaneously, three gases (CO₂, N₂, and O₂) were injected at 200 sccm into the bioreactor using the MCQ GB100 gas mixer (Monkey Industrial Supply, Irvine, CA, USA). The pH and dissolved oxygen (DO) in the bioreactors’ water were monitored and adjusted daily for the duration of the incubation to maintain a similar pH and DO to the data measured from the groundwater on the day of sampling. This was achieved by increasing or decreasing the amount of O₂ and CO₂ injected into the bioreactors. The whole bioreactor system was located in an incubation chamber, in which the temperature was set according to that measured in the groundwater during sampling (13 °C).

To test for external contamination during the experiments, a control run was carried out in the bioreactors. After cleaning the equipment, replacing the tubes, and sterilizing new rock chips, the assembled bioreactor systems were autoclaved. The incubation experiment was then run again in the same conditions detailed above, using sterile Milli-Q water instead of groundwater, for the same duration. Two rock chips from each control bioreactor were used for DNA extractions and sequencing, as explained below.

2.4. Subsample Collection during the Incubations in the Bioreactors and DNA Extractions

During the 24 days of incubation, 300 mL of excess water was filtered every two days from each of the three collection bottles to analyze the planktonic community using sterile polyethersulfone membrane filters with 0.2 μm diameter pores (Sartorius, Göttingen, Germany). Each obtained filter was stored at −20 °C. At the end of the experiment, all the rock pellets were taken from each of the 3 bioreactors, to study the colonization of rock surfaces (sessile community), and frozen at −20 °C. DNA was extracted from the rock surfaces of each rock pellet separately, as well as from the filters.

The DNA of the sessile microorganisms was extracted from each pellet following the protocols of Lazar et al. [14]. First, a 2x AE elution buffer was prepared (Supplementary Material Table S1). The rock chip was then placed in a 50 mL Falcon tube with 1 g of autoclaved 0.1 mm zircon/silica beads and 2 mL of 2x AE buffer, and vortexed for 30 s. Then, 500 μL of 20% SDS was added, and the solution was vortexed for 1 min at maximum speed. The solution was then centrifuged for 5 min at 8000× g, and the rock chip was put into a 15 mL Falcon tube. An amount of 2 mL of phenol/chloroform/isoamyl alcohol (25:24:1) was added, and the solution was vortexed at maximum speed for 1 min and then centrifuged for 5 min at 8000× g. The obtained aqueous phase was transferred to a new 50 mL tube, and 2 mL of chloroform was added. Then, the new tube was vortexed for 1 min and centrifuged for 10 min at 8000× g. The supernatant was then again transferred into a new 50 mL tube. To precipitate the DNA present in the sample, 2 μL of glycogen (20 mg/mL), 200 μL of sodium acetate (3M), and 4 mL of 100% ethanol were added to the aqueous phase. The tube was inverted and incubated overnight at 4 °C. The next day, the tube was centrifuged for 30 min at 8000× g. The obtained pellet was dried out and then rinsed with 0.5 mL of 70% ethanol. The tube was then centrifuged for 10 min at 8000× g, and the ethanol was completely removed. The pellet was then re-suspended in 100 μL of 1x Tris EDTA (TE) buffer.

DNA was extracted from the filters using the RNeasy PowerSoil Total RNA kit (QIAGEN, Hilden, Germany), followed by the Rneasy PowerSoil DNA elution kit (QIAGEN, Germany), following the instructions provided by the manufacturer. To test for contamination in the kit, we ran a blank extraction using MilliQ water filtered through the same filters used for the groundwater samples.

2.5. DNA Sequencing and Bioinformatics Analysis

Bacterial and archaeal 16S rRNA, and eukaryotic 18S rRNA genes from rock chips and filters were amplified separately using the polymerase Phusion Hot Start HF (ThermoFisher, Waltham, MA, USA). Primer pair B341F–B785R (CCTACGGGAGGCAGCAG and GACTACHVGGGTATCTAATCC) was used for Bacteria, A340F–A915R (CCCTACGGGGYGCASCAG and GTGCTCCCCCGCCAATTCCT) was used for Archaea, and E960F–E1438R (GGCTTAATTTGACTCAACRCG and GGGCATCACAGACCTGTTAT) was used for Eukaryotes (Supplementary Material Tables S2 and S3). Sequencing was carried out using the Illumina Miseq platform of the CERMO-FC (Centre d’Excellence en Recherche sur les Maladies Orphelines—Fondation Courtois) at UQAM, with the Miseq Reagent v3 600-cycle kit (Illumina, San Diego, CA, USA) and a paired reading of 300 bp. Negative PCR controls, as well as the control samples for water and rock samples, were also sequenced for all three domains. Despite all our efforts, we did not manage to amplify the archaeal 16S rRNA genes for the sessile community. Thus, these samples were not sequenced. The raw reads were deposited into the National Center for Biotechnology Information (NCBI) under the BioProject ID PRJNA1058723.

The DADA2 package in Rstudio v.4.2.2 was used to process the sequences and produce amplicon sequence variant (ASV) tables. Sequences present in the controls were considered contaminants and were removed from the datasets using the decontam package in R. For each domain (Bacteria, Archaea, and Eukaryotes), rarefaction was carried out using the median sequencing depth method. Taxonomy was assigned using the Silva 138.1 reference database.

2.6. Digital PCR Amplifications (dPCR)

To determine the absolute abundance of bacterial 16S rRNA genes on the rock surfaces, we performed dPCR amplifications of the sessile community. We used the QuantStudio 3D Digital PCR (ThermoFisher, Waltham, MA, USA) instrument as well as QuantStudio 3D PCR Master Mix v2 (ThermoFisher, USA) (Supplementary Material Tables S4 and S5). Results are expressed as gene copies/rock chip.

2.7. Statistical Analyses

Because the data was not distributed normally, non-parametric tests were used. The different water geochemical properties were compared between the bioreactors, as were the mineralogic characteristics and the bacterial sessile absolute abundances, using a Kruskal–Wallis test with the dunnTest function of the FSA package in R. The alpha diversity (Shannon index) was calculated using the diversity function from the vegan package (version 2.6.4). To compare the alpha diversity between both lifestyles (sessile and planktonic) for Eukaryotes and bacteria, as well as between the 3 bioreactors, we used a Kruskal–Wallis test.

Beta diversity was calculated using the Bray–Curtis dissimilarity matrix. We used PCoA ordination (principal coordinate analysis) with PAST v4 to observe the distribution and clustering of the different water and rock samples. We tested the influence of environmental parameters (bioreactor or time) on the β diversity indices of both the planktonic and sessile communities using PERMANOVA (permutational multivariate analysis of variance) with the adonis2 function of the vegan package in R. For the planktonic community, db-RDA (distance-based redundancy analysis) was performed to determine correlations between β diversity and the water geochemical properties with the capscale function of the vegan package in R. For the sessile community, db-RDA was performed to determine correlations between β diversity and the rock mineralogy. Variance partitioning of the significant environmental variables was assessed using the varpart function of the vegan package in R.

We identified shared and unique genera between groups of samples using the mothur software (v.1.47.0). For each domain, a linear discriminant analysis effect Size (LEfSe) was performed to identify which microbial taxa are primarily responsible for driving differences between previously determined significantly dissimilar groups. These analyses were performed using the mothur software with a one-against-one multiclass parameter when more than 2 sample groups were used. Finally, to determine the proportion of the planktonic community (sources) that contributed to the sessile community (sinks), we ran multiple fast expectation–maximization microbial source tracking (FEAST) [24].

3. Results

3.1. Geochemical Water Properties and Rock Characteristics

Overall, for all three bioreactors, four minerals were identified on the rock chip surfaces: mica, quartz, red (alkali) feldspar, and white (plagioclase) feldspar (Figure 3). There were also iron oxides present on the surfaces. Quartz was the most abundant mineral in all three bioreactors (19.41 to 98.14% of the surfaces).

In the in situ groundwater collected in the field, and subsequently used for the incubations, the pH was acidic (5.21), dissolved oxygen was at 54.5%, DIC and DOC were high (14.12 and 28.68 mg/L), and ammonia and nitrate concentrations were low but detectable (0.053 and 0.0045 mg/L).

DIC concentrations decreased drastically after 6 days of incubation, while DOC concentrations remained stable throughout the experiment, with peaks at 12, 14, and 18 days, depending on the bioreactor (Figure 4). Ammonia concentrations decreased during the experiment, with a peak after 10 days of incubation. Nitrate concentrations were low throughout the experiment, with an increase (peak) after 10 days of incubation and a second increase at the end of the incubations after 22 days. To check whether or not conditions within the bioreactors were similar, we compared the average surface mineral abundances, and all tested parameters were not statistically different between the bioreactors (Supplementary Material Table S6). The same was observed for all water geochemical variables (Supplementary Material Table S7).

3.2. Rock-Attached Bacterial Absolute Abundance and Correlation with the rock Properties

The bacterial absolute abundance numbers varied from 6.71 × 10⁴ to 2.39 × 10⁶ gene copies/rock chip for bioreactor 1, 1.94 × 10² to 6.53 × 10⁵ for bioreactor 2, and 1.43 × 10⁵ to 2.39 × 10⁶ for bioreactor 3 (Supplementary Material Table S8). Abundances were significantly different between bioreactors (BRs) 1 and 2, and between BRs 2 and 3 (p = 0.025 and 0.038, Kruskal–Wallis), but not between BRs 1 and 3 (p = 0.83, Kruskal–Wallis).

3.3. 16S and 18S rRNA Gene Taxonomic Composition

For the Bacteria domain, the peat bog acting as the primary water recharge source of the Covey Hill aquifer was dominated by Acidocella (10.33% of the total reads), Granulicella (7.95%), Legionella (7.38%), unclassified (unc.) Caulobacteraceae (6.92%), and unc. Acidimicrobiia (5.3%) (Figure 5a). For the Eukaryote domain, the peat bog was dominated by Crustacea (25.52%), Rhynchostomatia (19.94%), Peritrichia (10%), and Hypotrichia (5.1%) (Figure 5b). For the Archaea domain, the peat bog was dominated by Nanoarchaeia SCGC_AAA011-D5 (46.68%), Micrarchaeales CG1-02-32-21 (14.34%), Crenarchaeota group 1.1c (11.1%), and Nanoarchaeia GW2011_GWC1_47_15 (5.31%) (Figure 5c). The initial groundwater was composed of Sulfurimonas (35.59%), unc. Hydrogenophilaceae (6.35%), unc. Thermodesulfovibriona (6.091%), and unc. Acidimicrobiaa (5.86%) in terms of bacteria; Pezizomycotina (31.34%), Saccharomycotina (20.082%), and Hypotrichia (5.84%) in terms of Eukaryotes; and unc. Bathyarchaeia (54.66%), Nanoarchaeia SCGC_AAA011-D5 (19.94%), Crenarchaeota group 1.1c (7.063%), and Micrarchaeales CG1-02-32-21 (5.49%) in terms of Archaea.

For the bacterial planktonic community, we observed a similar temporal pattern in all three bioreactors (Figure 5a). Sulfurimonas was the most abundant genus in the beginning from days 2 to 10 and started decreasing in relative abundance after 2 days of incubation (Day 2, BR1: 51.74%; BR2: 28.37%; BR3: 30.1%). Pseudomonas and Collimonas were also important at the start of the experiment and started decreasing after 4 and 8 days of incubation. Meanwhile, Gaillonella become the dominant taxa after 12 days of incubation until 20 days (Days 12–20, BR1: 24.29–34.73%; BR2: 23.1–27.21%; BR3: 20.12–21.5%), followed by Alkanindiges and unc. Acetobacteraceae. For the sessile community, the rock samples were taken on the last day of the experiment, so no temporal variation could be analyzed. For all samples from all three bioreactors, the same three genera dominated: Gaillonella (BR1: average of 25.73%; BR2: 19.49%; BR3: 28.55%), Alkanindiges (BR1: 14.44%; BR2: 19.4%; BR3: 16.65%), and unc. Acetobacteraceae (BR1: 13%; BR2: 12.47%; BR3: 9.75%). The other identified taxa (representing more than 2% of the total community) were Collimonas, Sulfurimonas, Sideroxydans, unc. Hydrogenophilaceae, unc. Acidimicrobiia, Rhodoferax, Afipia, Undibacterium, Methylotenera, Rhodoblastus, Cupriavidus, Herminiimonas, Rhodovastum, unc. Oxalobacteraceae, and Acidocella.

In the eukaryotic planktonic community, we observed a temporal pattern in all three bioreactors as well (Figure 5b). Pezizomycotina was the most abundant genus after 2 days of incubation (Day 2, BR1: 62.76%, BR2: 78.75%, BR3: 45.32%), but rapidly decreased in relative abundance after that. Unc. Chrysophyceae was the second most dominant taxa after 2 days of incubation, and dominated after 4 days (Day 4, BR2: 25.4%, BR3: 47.32%). After 6 days, Synurales dominated the planktonic community, until 12 days into the experiment. Finally, after 14 days for BR2 and after 18 days for BR1 and BR3, Apoikiales dominated the planktonic community. Synurales, Eimeriida, and Glissomonadida were also an important part of the community at the later time points. For the sessile community in bioreactors 2 and 3, the same three genera dominated: Apoikiales (BR2: average of 37.9%; BR3: 21.66%), Glissomonadida (BR2: 27.68%; BR3: 40.91%), and Synurales (BR2: 16.36%; BR3: 17.39%). Bioreactor 1 was dominated by two of these taxa: Apoikiales (34.65%) and Synurales (31.11%). The other identified taxa (representing more than 2% of the total community) were unc. Chrysophyceae, Pezizomycotina, Choanoflagellatea, Peritrichia, Rozellomycota, and unc. Amoebozoa

For the Archaea domain, due to the inability to sequence any of the DNA extracted from the rock chip surfaces, we only analyzed the planktonic community (Figure 5c). Unlike the two other domains, we did not observe major temporal patterns at the genus level in the planktonic community. Most of the dominant taxa identified in the initial groundwater sample were observed in the planktonic samples (unc. Bathyarchaeia, Nanoarchaeia SCGC_AAA011-D5, and Micrarchaeales CG1-02-32-21, as well as Nanoarchaeia GW2011_GWC1_47_15). We also detected the following genera: Crenarchaeota group 1.1c, candidatus (cand.) Micrarchaeum, unc. Methanomassiliicoccales, and cand. Nitrosotalea.

For the Bacteria, overall, the sessile and planktonic communities shared 57.05% of the genera; 16.67% were unique to the planktonic community, and 26.28% were unique to the sessile community (Supplementary Table S9). For Eukaryotes, the sessile and planktonic communities shared 56.25% of the genera, with 35.71% being unique to the planktonic community, and 8.04% being unique to the sessile community.

3.4. Alpha Diversity Indices

For the Bacteria domain, the α diversity indices of the sessile community varied from 2.334 to 3.864, and from 2.894 to 4.081 for the planktonic community (Supplementary Figure S2). There was a significant difference in the diversity indices between both lifestyles, with indices for the sessile community being higher on average (Supplementary Table S10). There was a significant difference in the sessile α diversity indices between bioreactors 1 and 3. There was no significant difference in the planktonic α diversity indices between the 3 bioreactors. For the Eukaryote domain, the α diversity indices varied from 1.911 to 3.079 for the sessile community, and 1.884 to 3.636 for the planktonic community (Supplementary Figure S3). There was no significant difference in the diversity indices between both lifestyles, or between bioreactors for both the planktonic and sessile communities (Supplementary Table S10). For the Archaea domain, the planktonic α diversity indices varied from 4.683 to 5.272 (Supplementary Figure S4). There was no significant difference in the diversity indices between bioreactors (Supplementary Table S10).

3.5. Beta Diversity Ordination Analyses and Influence of Environmental Variables

A PCoA was run for each domain to view the β diversity-based similarities between samples. For this analysis, we included the initial in situ groundwater and peat bog samples collected in the field on the day of sampling. For the Bacteria domain, the PCoA plot showed that the sessile community samples were all grouped based on which bioreactor they were incubated in, especially for samples incubated in bioreactor 1 (Figure 6a). The sessile community samples were most similar to the planktonic samples from the later incubation days (14 to 22). The planktonic community samples were grouped based on incubation time: day 2 samples clustered together, as did samples from days 4 to 12 and samples from days 14 to 22. The in situ groundwater samples were closest to the water from day 2. To support these observations, we ran PERMANOVA analyses, all excluding the field samples. The first one included both the sessile and planktonic communities, indicating that lifestyle was significantly correlated with the bacterial community and explained 29.95% of the variance (

Table 1. PERMANOVA analyses of the sessile and planktonic communities from the three domains. SumOfSqs, sum of squares.

Domain	Clusters	Parameters	df	SumOfSqs	R2	F	p-Value
Bacteria	Lifestyle (sessile/planktonic)	Type	1	2.5533	0.29953	26.939	0.001
		Residual	63	5.9710	0.70047
		Total	64	8.5243	1.00000
	Sessile	Bioreactor	2	0.72573	0.46152	13.713	0.001
		Residual	32	0.84675	0.53848
		Total	34	1.57248	1.00000
	Planktonic	Time group	2	2.5023	0.56888	21.2485	0.001
		Bioreactor	2	0.4243	0.09646	3.6027	0.008
		Residuals	25	1.4720	0.33466
		Total	29	4.3985	1.00000
Eukaryote	Lifestyle (sessile/planktonic)	Type	1	1.2931	0.17545	8.5111	0.001
		Residual	40	6.0771	0.82455
		Total	41	7.3702	1.00000
	Sessile	Bioreactor	2	0.69439	0.38606	5.345	0.004
		Residual	17	1.10427	0.61394
		Total	19	1.79867	1.00000
	Planktonic	Time group	2	2.5194	0.58885	14.0797	0.001
		Bioreactor	2	0.2382	0.05566	1.3309	0.219
		Residuals	17	1.5210	0.35549
		Total	21	4.2785	1.00000
Archaea	Planktonic	Time group	2	0.46131	0.28350	3.3878	0.001
		Bioreactor	2	0.21270	0.13072	1.5620	0.021
		Residuals	14	0.95318	0.58578
		Total	18	1.62719	1.00000

). A second PERMANOVA was carried out on the sessile communities, using the bioreactor as an environmental variable. Here, the bioreactor was a significant parameter, explaining 46.15% of the variance (Table 1). A third analysis was carried out on the planktonic communities using the group of incubation days identified in the PCoA plot and the bioreactor as environmental variables. Both variables were significantly correlated with the bacterial community, explaining 56.69% of the community variance for the incubation time, and 9.65% of the variance for the bioreactor (Table 1).

For the Eukaryote domain, the PCoA plot showed that the sessile community (Figure 6b), as observed for the Bacteria, was most similar to the planktonic community at later incubation time points (days 14–22). The planktonic community samples could be clustered into three groups based on incubation time: 2–4, 6–12, and 14–22 days. The in situ groundwater samples were closer to the sample incubated for 2–4 days. The PERMANOVA analysis indicated that that lifestyle (sessile/planktonic) was significantly correlated with the eukaryotic community and explained 17.54% of the variance (Table 1). For the sessile community, the bioreactor was a significant parameter, explaining 38.61% of the variance. For the planktonic communities, only the incubation time parameter (and not the bioreactor) was statistically significant, explaining 58.89% of the community variance.

For the Archaea domain, the PCoA plot showed one cluster containing most of the field in situ groundwater samples and the planktonic community samples from 2–8 days of incubation (Figure 6c), while a second cluster contained samples from days 10 to 16, and a third cluster contained samples from days 18 to 20. The PERMANOVA analysis indicated that both the incubation time and bioreactor were significantly correlated with the archaeal community (Table 1), explaining 28.35% and 13.1% of the variation.

3.6. Community Structure Correlation with the Water Geochemical and Rock Properties

We ran db-RDA analyses using the same β diversity datasets as those used for the PERMANOVA and excluding the in situ groundwater and peat bog water samples. The planktonic dataset was analyzed separately from the sessile dataset. For the sessile community, the rock mineral properties were used as environmental variables, while for the planktonic community, the water geochemical parameters were used. For the bacterial sessile community, we added bacterial absolute abundances (measured with dPCR) as an environmental variable, and the db-RDA showed that quartz, white feldspar, and bacterial absolute abundances were significantly correlated with beta diversity (Supplementary Material Table S11). The graph indicated that quartz was correlated with bioreactor 1, while white feldspar was correlated with bioreactors 2 and 3 (Figure 7a). Variance partitioning showed that quartz explained 1.28% of the community variance, white feldspar explained 0.37% of it, and bacterial abundance explained 0.2%. Furthermore, the db-RDA showed that only DIC was significantly correlated with beta diversity of the planktonic community (Supplementary Material Table S11), and the graph indicated that DIC was correlated with later time point samples (18–22 days) (Figure 7b). Variance partitioning showed that DIC explained 6.24% of the community variance.

For the Eukaryote domain, none of the mineral variables were significantly correlated with the beta diversity of the sessile community (Supplementary Material Table S11). The db-RDA also showed that only DIC was significantly correlated with the beta diversity of the planktonic community (Supplementary Material Table S11), and the graph indicated that DIC was correlated with later time point samples (18–22 days) (Figure 8). Variance partitioning showed that DIC explained 12.94% of the community variance.

For the Archaea domain, the db-RDA also showed that both DIC and DOC were significantly correlated with the beta diversity of the planktonic community (Supplementary Material Table S11), and the graph indicated that both variables were correlated with later time point samples (18–20 days; DOC for bioreactor 2, and DIC for bioreactors 1 and 3) (Figure 9). Variance partitioning showed that DIC explained 9.34% of the community variance, and DOC explained 10.53% of it.

3.7. Discriminative Genera Explaining Differences between Community Groups

In the LEfSe analyses, we looked at which genera for both the Bacteria and Eukaryote domains were significantly different in terms of relative abundance for both lifestyles (planktonic and sessile). For the sessile bacterial community, the genera with the highest LDA (linear discriminant analysis) scores were Gallionella, Alkanindiges, Rhodoferax, Undibacterium, Methylotenera, and Sideroxydans, (Supplementary Material Table S12). For the planktonic community, the genera with the highest LDA scores were Sulfurimonas, Collimonas, and Pseudomonas. For the sessile eukaryotic community, the genera with the highest LDA scores were Glissomonadida, Apoikiales, unc. Amoebozoa, Peritrichia, Rhynchostomatia, and unc. Alveolata (Supplementary Material Table S12). For the planktonic community, the genera with the highest LDA scores were Eimeriida, Saccharomycotina, Choanoflagellatea, and Rozellomycota. Furthermore, we used LEfSe to determine discriminative genera in the different bioreactors for the bacterial and eukaryotic sessile communities. For the Bacteria domain, for BR1, the genera with the highest LDA scores were unc. Hydrogenophilaceae, Rhodoferax, unc. Burkholderiaceae, and unc. Paludibacteraceae; for BR2 these were Undibacterium, Methylotenera, Occallatibacter, Collimonas, Cupriavidus, and Variovorax; and for BR3, these were Rhodoblastus, Sulfurimonas, Novosphingobium, and Acidocella (Supplementary Material Table S13). For the Eukaryote domain, for BR1, the genera with the highest LDA scores were Synurales and Peritrichia; for BR2 and BR3, there was no significant genus (Supplementary Material Table S13).

Finally, we used LEfSe to distinguish the discriminative genera in the different incubation time groups for all three domain planktonic communities. For the Bacteria domain, for day 2, the genera with the highest LDA scores were Sulfurimonas, and unc. Acidimicrobiia; for days 4–12, these were Collimonas, Undibacterium, and Sideroxydans; and for days 14–22, these were Gallionella, unc. Acetobacteraceae, and Alkanindiges (Supplementary Material Table S14). For the Eukaryote domain, for days 2–4, the genera with the highest LDA scores were Pezizomycotina, Saccharomycotina, and Embryophyceae; for days 6–12, Synurales had the highest LDA score; and for days 14–22, Apoikiales had the highest score (Supplementary Material Table S14). For the Archaea domain, for days 2–8, the genus with the highest LDA score was unc. Bathyarchaeia; for days 10–16, there was no significant genus; and for days 18–20, SCGC AAA011_D5 Woesearchaeales, GW2011_GWC1_47_15 Nanoarchaeia, unc. Woesearchaeales, and CG1_02_32_21 Micrarchaeales were the genera with the highest LDA scores (Supplementary Material Table S14).

3.8. Microbial Source Tracking

Microbial source tracking analyses were carried out for the Bacteria and Eukaryote domains. For the rock-attached sessile communities, we used the peat bog, the in situ groundwater, and the water samples from the previous time points as potential sources. For the bacterial sessile community, for which we only have the final time point, the water samples from days 14 and 16 were the biggest contributors (Figure 10a). The water samples from the 22nd day of the experiment (the latest water timepoint) did not contribute as much to the sessile community as the time points mentioned previously. The in situ groundwater contributed, on average, 2.29% for bioreactor 1, 0.68% for bioreactor 2, and 1.94% for bioreactor 3. For the eukaryotic sessile community, the water samples from days 14, 16, and 18 were the biggest contributors (Figure 10b). As observed for the Bacteria, the water samples from the 22nd day did not contribute as much to the sessile community as the time points mentioned previously. The in situ groundwater contributed, on average, 2% for bioreactor 1, 0.24% for bioreactor 2, and 3.18% for bioreactor 3.

4. Discussion

4.1. Evolution of Planktonic Communities over Time

Incubation time had a significant effect on the variation in the bacterial, eukaryotic, and archaeal planktonic composition, explaining a higher variance in the community composition compared with the bioreactor factor. This shows that the microbial community structure changed during the 24 days of the experiment, probably as a result of the water geochemical parameter changes and rock surface colonization of some of these planktonic microorganisms. Hence by leaving the planktonic community to settle on the surfaces, they most likely changed the community composition and structure.

For the Bacteria domain, Sulfurimonas was significantly higher in relative abundance in the water at the beginning of the experiment (first 2 days) and was also the dominant taxon in the initial groundwater used for the incubations, but its relative abundance decreased after 4 days. Sulfurimonas is a genus mostly found in sulfidic habitats such as hydrothermal vents, marine sediments, or terrestrial oil fields [25]. These chemolithoautotrophic bacteria use sulfide, elemental sulfur, and thiosulfate to generate energy, and CO₂ as a carbon source [26], which was measured in high amounts in the initial groundwater. Although we did not measure sulfur compounds in the water, there probably was sulfide in the in situ groundwater based on the smell of sulfur that could be detected while sampling; therefore, this bacterium would have been able to thrive. One reason that might explain the decline in this genus in the experiment could be a reduction in sulfur or any sulfur derivative in the water during the incubations, as well as the observed sharp DIC decrease after 6 days.

Then, Collimonas, Undibacterium, and Sideroxydans were significantly higher in relative abundance in water samples from days 4 to 12. Collimonas is a chitinolytic soil bacterium isolated from slightly acidic dune soils and is able to grow on living fungal hyphae [27]. Whole-genomic sequencing suggested its capacity for swimming motility and siderophore production (iron acquisition when it is in limited amounts) [28]. Undibacterium was isolated from drinking water [29], and U. oligocarboniphilum was shown to thrive in low-carbon-substrate concentrations [30]. This genus also has the ability to oxidize iron, as demonstrated in groundwater-fed sand filters [31]. Finally, Sideroxydans is a chemolithoautotrophic iron oxidizer, found in freshwater, and often in acid habitats [32,33]. Thus, the major bacterial taxa found in the water during the midpoint of the experiment were likely coping with oligotrophic conditions, and were all able to oxidize iron, which likely came from the rock chips, since we detected iron oxides on all rock surfaces, and since it is also found in mica and feldspar.

Finally, the later time points (days 14 to 22) were conducive to Gallionella, Alkanindiges, and unc. Acetobacteraceae, which were present in significantly higher amounts at these time points compared with the other incubation times. Gallionella is an acid-tolerant iron-oxidizing lithoautotrophic bacterium found in soils [34], groundwater [35], and low-nutrient aquatic environments rich in calcium [36]. Alkanindiges is an obligate hydrocarbonoclastic bacterium that was first isolated from oilfield soils [37] and was also detected on a marble statue [38]. Acetobacteraceae is an important family of bacteria, used for industrial fermentation in the food sector. It is composed of two major groups: acetous and acidophilic bacteria. The acetous group comprises mainly acetic acid-fermenting bacteria. Acetobacteraceae are important secondary metabolite producers [39,40], and can also fix nitrogen [41]. The activity of these bacteria probably decreased the pH in the water, leading to the establishment of acid-tolerant iron-oxidizers like Gallionella. The increase in DIC, significantly correlated with the planktonic community structure in the later incubation days, also explains the dominance of this lithoautotrophic genus. The acidophilic family Acetobacteraceae contains genera like Acidiphilium that grow in organic-poor media between pH 1.9 and 5.9 [42]. If the more labile carbon sources were first used during the first days of the experiment, this would explain the detection of bacteria able to use more complex carbon sources after 14 days of incubation.

For the eukaryotic planktonic community, Pezizomycotina, Saccharomycotina, and Embryophyceae were significantly higher in relative abundance during the first 4 days of the experiment and were also major taxa in the in situ groundwater. Pezizomycotina are fungi of the Ascomycota phylum [43]. They are a highly diverse group including saprophytes, plants, and mutualists. They have been found in aquifer ecosystems [44], where they are surmised to play a role in nutrient cycling. The fungal subphylum Saccharomycotina is made up of yeasts with very diverse nutritional metabolisms inhabiting major aquatic and terrestrial habitats [45]. The Embryophyceae were likely an initial food source for these fungi, and predation by the bacterial genus Collimonas might explain their decline after 4 days.

From 6 to 12 days, Synurales were significantly higher in relative abundance in the water samples. These Eukaryotes are unicellular, motile, siliceous scale-bearing heterokont algae found mostly in freshwater ponds and lakes [46]. After 14 days, Apoikiales dominated the water. Apoikia lindahlii was isolated from acidic lakes, and these are biflagellate bacterial predators; they also consume dissolved organic molecules [47]. As such, they are considered important players in aquatic food webs. The Apoikiospumella are considered Spumella-like flagellates [48]. A study conducted by Rothhaupt [49] looked into the competition for nutrients between Ochromonas algae and Spumella, showing that Spumella were better at acquiring food than Ochromonas in a dark environment. The authors mentioned that Ochromonas have difficulty growing in the dark and need a sufficient concentration of bacteria to feed on. A similar competition outcome between Synurales and Apoikiales could explain the temporal shift that we observed in both groups over time. As for the Bacteria domain, the dominance of (autotrophic) algae at the later time points can be linked to the increase in DIC, which was significantly correlated with the planktonic community structure. It is highly probable this DIC increase was the result of metabolic activities carried out by the fermenting and heterotrophic bacteria and fungi during the first phases of the bioreactor experiment.

For Archaea, we identified three temporal groups in the water during the 24 days of the experiment. From 2 to 8 days, Bathyarchaeia was the dominant lineage, and this group contains some of the most ubiquitous and abundant Archaea on Earth, with very diverse metabolisms, such as in those for C₁ compound or lignin utilization [50], all of which are probable molecules seeping down from the surface peat bog water recharge. The later time points (18–22 days) were dominated by uncultured lineages. SCGC AAA011_D5 Nanoarchaeia (Woesearchaeales) were found in subarctic rivers, where they were linked to high DOC concentrations [51], and in groundwater [52], and were involved in kelp degradation after 60 days when most of the labile organic carbon had been used up [53]. Genomic reconstructions suggest that the SCGC_AAA011-D5 Archaea are organic compound scavengers. GW2011_GWC1_47_15 Nanoarchaeia have been found in sludge [54] and dump site soils [55]. Nanoarchaeia seem to be linked to the degradation of molecules with higher DOC aromaticity [56]. CG1_02_32_21 Micrarchaeales played a role in the co-compositing of food and agricultural waste [57]. These observations point to complex organic carbon utilization and processing by Archaea during the last days of incubation, as supported by the significant correlation with DOC concentrations.

4.2. Microbial Colonization of the Rock Surfaces

After 24 days of incubation in the bioreactors, more than half of the genera were shared between the planktonic and sessile communities for the Bacteria and the Eukaryotes. Since the rock surfaces were sterile at the start of the experiment, this suggests that many aquifer microorganisms can switch lifestyles from planktonic to sessile, but probably also from sessile to planktonic [58]. Nonetheless, our study showed that bacterial α diversity indices were significantly different between both lifestyles, with sessile indices being higher, while eukaryotic diversity indices did not significantly differ. This could mean that the sessile lifestyle is preferred by many bacterial genera over the planktonic lifestyle. Indeed, colonization and development on rock surfaces is quite common in aquifers [59,60], and the colonization of minerals on surfaces can offer physical protection against the water flow [61]. The β diversity indices were also significantly different between both lifestyles for both the Bacteria and Eukaryote domains, explaining more than a quarter of the community variance for the Bacteria.

We found that the genera that were significantly higher in relative abundance in the bacterial sessile community potentially have different metabolisms. Some genera are involved in the iron cycle: Gallionella and Sideroxydans are iron-oxidizers, and a species within the Rhodoferax genus is an iron reducer [62]. Undibacterium was isolated from drinking water [29], and one species, U. oligocarboniphilum, was shown to thrive under low-carbon conditions [30]. Mehylotenera is an obligate methylamine-utilizing bacterium [63], and also uses methylated amines [64]. All taxa that were significantly higher in relative abundance in the sessile community were shown to be involved in biofilm communities. Although we did not check whether or not our sessile microorganisms did indeed form biofilms on the rock surfaces, another recent study showed, using SEM, that aquifer microbes form biofilms on solid surfaces [18]. Alkanindiges was involved in biofilm formation on metal surfaces [65], Undibacterium and Rhodoferax were prevalent in biofilms formed on decaying leaves [66], and Undibacterium was surmised to participate in secondary colonization of biofilm in a drinking water distribution system, as supported by its capability to use secondary metabolites produced by early colonizers [67]. Undibacterium and Methylotenera can form biofilms on microplastics that are frequently present in wastewater [68,69,70]. Methylotenera can participate in denitrifying biological filter biofilms emended with methanol [71]. Finally, Rhodoferax participated in riverbed biofilm formation after a sewage fungus outbreak [72]. Hence, numerous potential bacterial heterotrophic metabolisms seem to occur within the biofilm, exhibited by microbes adapted to living in low-carbon environments, or using C₁ compounds as carbon sources.

For the Eukaryote domain, the cercozoan biflagellate bacterivorous genus Glissomonadida was significantly higher in relative abundance in the sessile community. These protozoa are mostly found in soils [73], where they graze on bacterial biofilms [74]. In fact, the formation of biofilms is thought to be a bacterial strategy used to avoid predation by protozoa [75]. Peritrichia were also significantly higher in relative abundance in the sessile community. These are ciliated protists living in aquatic habitats and soils [76]. A study looking at biofilm formation during wastewater treatment showed that Peritrichia settle on sludge flocs and build stalks, which are subsequently used and colonized by bacteria that end up killing the protists [77]. Therefore, it is unclear whether or not these Eukaryotes are actually part of biofilm functioning, or if they are associated with the biofilm through its formation and predation.

Microbial source tracking for the sessile communities in both the bacteria and Eukaryote datasets showed that the dominant sources were the planktonic samples from the later days of the experiment (14 to 18 days) but not the last days. This suggests that the colonization of rock surfaces by planktonic microbes occurred up until late into the experiment, but also that after 18 days, the exchange between both communities decreased. Competition among microorganisms for resources, hiding from Eukaryote predation, and space on the surfaces [78] could explain why organisms from the planktonic community continued to contribute to the sessile community until the final week of the incubations. Microbial primary succession, during which early pioneers could have then been replaced by other taxa [79,80], might have led to the stabilization of the sessile community in the final week of the experiment.

For the Archaea domain, it was not possible to retrieve sequences from the rock surface samples. Archaeal species have previously been detected in biofilm samples, but the environmental conditions of the aquifers were different from those in the present study [81]. The low pH in the water, the characteristics of the rock surfaces, and the presence of eukaryotic grazers could all potentially explain the absence of Archaea in the sessile community. It is also possible that the Archaea initially present in the groundwater did not possess the necessary biochemical pathways to attach to solid surfaces, carry out quorum sensing, or produce extracellular polymeric substances [82,83].

4.3. Influence of the Rock Properties on Surface-Attached Bacterial Communities

Quartz, white feldspar, and bacterial absolute abundances were significantly correlated with the bacterial sessile community structure, although they explained less than 3% of the community variance. Previous studies have found that some bacteria could colonize quartz [84,85] and another study conducted on soil found that feldspar (more precisely potassium feldspar) also influenced the microbial community [86]. It is possible that these bacteria can utilize some compounds found in minerals for their growth and survival, such as phosphorus, calcium, and magnesium, explaining their prevalence on the rock surfaces. Furthermore, the bioreactor in which the incubations were carried out also had a significant influence on the sessile bacterial community structure, explaining almost half of the community variance. The average surface mineral abundances were not significantly different between the three different bioreactors. Thus, the differences in community composition between the bioreactors could be explained by local differences in the water’s geochemical composition, since the rock surfaces and communities were continuously exposed to the water flowing in the system. Seepage of groundwater and microbial activities—the occurrence of which was supported by differences in the bacterial absolute abundances between bioreactors—can modify the pH and water geochemistry, and could have altered the mineral composition of the rocks over the course of the experiment. Local differences for each rock chip in each bioreactor could have contributed to the observed influence of the bioreactor on the sessile community composition.

The LEfSe analysis indicated that specific bacterial genera dominated the sessile communities from each bioreactor. Rhodoferax was higher in relative abundance in bioreactor 1, and had a community composition correlated with that of quartz. This putative hydrogenotrophic iron-reducing taxon was shown to be influenced by the presence of iron minerals in a basaltic glacier in Iceland [87]. Unidbacterium, Collimonas, and Cupriavidus were higher in relative abundance in bioreactor 2, while Acidocella was higher in relative abundance in bioreactor 3, and had community compositions correlated with those of white feldspar. Undibacterium was part of the community involved in regolith-hosted rare-earth-element deposits [88], as well as the biomineralization of stone surfaces containing calcium carbonates and phosphates [89]. Basalt leaching of nutrients such as calcium, iron, phosphorus, or magnesium supported the growth of Cupriavidus metallidurans [90]. Collimonas pratensis is very efficient in weathering minerals such as biotite, olivine, garnet, hematite, and apatite, and its metabolic capacities are linked to the buffering capacity of the environmental conditions [91]. Finally, Acidocella was linked to the mineralogical composition (pyrite and quartz) in an acidic metal mine [92]. Thus, it is possible that the local heterogeneity of the mineral composition in each bioreactor influenced colonization by the different taxa depending on metabolic capabilities and uptake efficiency. Moreover, Jones et al. demonstrated that quartz was more stable in acidic conditions, while feldspar mobilized potentially toxic metals like aluminum [11]. They showed that neutrophilic bacteria preferentially colonize high-buffering carbonate rocks, while acidophiles prefer non-buffering minerals like quartz. It is hence likely that the pH conditions found in the water surrounding the rock surfaces also impacted which bacteria attached the surfaces, and on which mineral types.

5. Conclusions

In summary, we observed a temporal variation among planktonic communities for all studied domains (Bacteria, Eukaryotes, and Archaea). There was a significant difference in α diversity indices between lifestyles (sessile/planktonic) for bacterial communities, but not in Eukaryotes. The mineral composition of the rock surfaces had a significant influence on the β diversity for the Bacteria domain, with quartz and feldspar influencing the community structure. Many bacteria involved in iron cycling were found in both the water and rock surfaces but did not dominate the initial in situ groundwater, suggesting that rock colonization and subsequent sessile microbial activities enhanced iron scavenging and weathering from the minerals. Our data also indicated probable predator-prey interactions between Bacteria and Eukaryotes.