For microbial analyses, about 20 L of water was passed through a 0.2 μm membrane filter paper (Pall Gelman Laboratory, Ann Arbor, MI) using stainless filter holder (Millipore Corp.) for 30 min, which was connected to the peristaltic pump in the springs or to the discharge ports in the wells. After filtration the filter papers were stored in a 50 mL conical tube and then stored at −20°C in the lab. All equipment was sterilized in an autoclave before sampling.

Samples were frozen at −20°C until the experiments were completed. The filters were cut into several small pieces and dispensed into a sterile microcentrifuge tube. Total genomic DNA was extracted by using an i-geonomic Soil DNA Extraction Mini Kit (iNtRON, South Korea) with a bead-beating apparatus, according to the manufacturer's directions. DNA concentration was quantified with a Qubit fluorometer (Invitrogen, USA), by following the manufacturer's instructions. The extracted DNA was used for PCR amplification for clone library and MiSeq sequencing as described below.

For PCR amplification, the DNA template (100 ng) was added to a master mix of ddH₂O, 2.5 μL of 10 × PCR buffer, 1 μL of 50 mM MgSO₄, 0.5 μL of 10 mM dNTPs, 2 μL of 8F primer (5′-AGAGTTTGATCMTGGCTCAG-3′)(10 μM), 2 μL of 1492R primer (5′-TACGGYTACCTTGTTACGACTT-3′)(10 μM), and 0.2 μL of i-Taq DNA polymerase (5 U μL⁻¹, Intron, South Korea), and brought to a final volume of 25 μL with sterile ddH₂O. For amplification of archaeal DNA, primers 109F (5′-ACKGCTCAGTAACACGT-3′) and 958R (5′-YCCGGCGTTGAMTCCAATT-3′) from the same source were used. Templates were amplified by using a thermal cycler. The initial denaturation step at 94°C for 2 min was followed by 30 cycles of 94°C (30 s), 58°C (45 s), and 72°C (90 s), with a final extension at 72°C for 10 min. The PCR products were purified as indicated above. Clone libraries were constructed by using the T-Blunt™ PCR Cloning Kit (with DH5α Competent Escherichia coli; Solgent, South Korea). Approximately 50 colonies per plate were selected randomly and re-amplified by using the vector-specific primers (0.25 μM each) M13F (5′-GTAAAACGACGGCCAG-3′) and M13R (5′-CAGGAAACAGCTATGAC-3′) (10 μM, BIONEER, Korea). M13 primer-amplified PCR products were visualized on a 1% agarose gel (Invitrogen) to check for successful amplification of the cloned inserts. The successful products were then treated with ExoSAP-IT (USB, CA, USA) (1.5 μL of product, 3.25 μL of ddH2O, and 0.25 μL of ExoSAP-IT for each reaction). Purified M13 PCR products were sequenced at the Marcrogen (Seoul, Korea) with the 8F forward primer.

All sequences were trimmed by using DNA Baser (Heracle BioSoft, Romania) to remove inconclusive nucleotides from the ends, on the basis of the number of ambiguous base calls and quality scores. The sequences were also screened for chimeras and other artifacts with the Mallard program (http://www.cardiff.ac.uk/biosi/research/biosoft/) (Ashelford et al., 2006). A multiple-sequence alignment prepared for each included the sequence E. coli U00096 as a reference. The sequences screened were analyzed vs. publicly available bacterial 16S rRNA gene sequences from the Ribosomal Database Project (RDP) Release 11.5 consists of 3,356,809 aligned and annotated 16S rRNA sequences through the RDP-II website (http://rdp.cme.msu.edu/). Distance matrices based on alignment of the library sequences to the RDP database were used as input for further analysis with the Distance Based Operational Taxonomic Units (OTUs) and Richness Determination (DOTUR) program (Schloss and Handelsman, 2005). All clones having a mutual sequence similarity of more than 97% were grouped into OTUs by using DOTUR. Rarefaction curves (number of OTUs detected vs. the number of clones analyzed) were calculated to explore the diversity of microbial communities.

We used two sets of primer pairs for MiSeq library construction. Taxonomic coverage of two forward primers and two reverse primers was evaluated by probe match tool of the Ribosomal Database Project (RDP) according to previous studies (Cole et al., 2014; Hugerth et al., 2014). For the first library, we used the primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 805R′, (5′-TACHVGGGTATCTAATCC-3′) (Caporaso et al., 2011) targeting hypervariable region 4 (V4) of the 16S rRNA genes of bacteria and archaea. The PCR conditions were 94°C for 3 min, followed by 30 cycles of 94°C for 30 s, 50°C for 60 s, and 72°C for 90 s, followed by a final extension at 72°C for 10 min. For the second library construction, the DNA extracted from groundwater samples were amplified with 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′) set of primers (Klindworth et al., 2012) targeting the V3–V4 of 16S rRNA gene. The PCR conditions were 95°C for 5 min, followed by 30 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s, followed by a final extension at 72°C for 5 min. Each sample was amplified in three replicate 25-μL PCR reactions and then purified using the UltraClean PCR clean-up kit (Mo Bio Laboratories, Solana Beach, CA, USA) and combined into a single volume. Purified DNAs were quantified by incorporating Picogreen (Invitrogen) according to the manufacturer's instructions. 16S rRNA-based MiSeq sequencing with a paired-end read length of 2 × 300 bp was conducted at Macrogen (Seoul, Republic of Korea) using the Illumina MiSeq system (Illumina, USA), according to the manufacturer's instructions.

To perform quality filter and merge simultaneously, we used the Illumina Utilities library as described previously in Eren et al. (2013) with an enforced Q30-check and ≤3 mismatches required to retain a sequence and its source code is available from (https://github.com/meren/illumina-utils). Sequencing data were analyzed using QIIME 1.8.0 (Caporaso et al., 2010). The detailed QIIME parameters used in the current analysis are following: The “pick_otus.py” script was used to cluster reads from all samples into OTUs based on nucleotide sequence identities. The clustering method was “Uclust” (Edgar, 2010). As default parameters, the similarity threshold was set at 0.97 (roughly corresponding to species-level OTUs). After picking OTUs, the “pick_rep_set.py” script was used for picking a representative set of sequences. For the best taxonomic assignment of the sequences, we used two taxonomy reference databases, RDP and SILVA. Both databases showed similar results, but SILVA database showed better taxonomic assignment particularly for Archaea (Table S1). The “assign_taxonomy.py” was used to assign the taxonomy of each sequence, with the 24 July 2014 release of the ARB-formatted SILVA small subunit (SSU) reference database (SSU Ref v.119) downloaded from www.arb-silva.de/documentation/release-119 (Quast et al., 2012). The “align_seqs.py” was used to align with the SILVA 119 database, and the template file used was SILVA 119 reference “core_Silva119_alignment.fna.” To filter out chimeric sequences after the alignment “filter_alignment.py” script was used with the template (lanemask_in_1s_and_0s). The “make_phylogeny.py” script was used with the “FastTree” program (Price et al., 2010) to construct a phylogenetic tree. The “make_otu_table.py” script was used to generate an OTU table in the biom format. The “filer_otu_from_otu_table.py” script was used to remove all OTUs that were observed fewer than two times (i.e., singletons). The “filter_taxa_from_otu_table.py” script was used to remove Archaeal taxa from an OTU table.

To calculate species diversity and richness within individual samples, alpha diversity analysis, including observed OTUs, Chao1 (Chao, 1984), Shannon (Shannon, 2001), Simpson (Simpson, 1949), ACE (Chazdon et al., 1998) indices for statistical tests, was processed using the QIIME script. Good's coverage estimator was also used to calculate the sequence coverage obtained for the 16S rRNA datasets (Good, 1953). To measure the similarity among communities, beta diversity was calculated based on weighted and unweighted unifrac distance matrices between communities and then 3D principle coordinates analysis (PCoA) plots were constructed. Non-metric multidimensional scaling (NMDS) analyses were performed with the function “meta.mds()” of the package vegan implemented in R (R Core Team, 2014), based on the microbial community information (relative abundance of each OTU detected by either MiSeq sequencing or clone library). The vector fitting of environmental variables to the NMDS ordination was determined by the “envfit()” function in the vegan package with 16 major components of physical and chemical characteristics. Significance was determined based on Bray–Curtis distances and 10,000 random permutations. Diversity metrics were calculated based on the sequencing data analyzed by 341-805R, with the minimum number of reads being 30,714 for bacteria and 1,974 for Archaea (Tables S1, S2).

Sequencing reads produced from Miseq during this study have been deposited at the NCBI Short Read Archive(SRA) under project accession number SRP119027 (Accession number: SRR6110859-SRR6110865).

Combining microbial community investigation with groundwater chemistry in subsurface environments is essential to better understand our ecosystem because many subsurface biogeochemical processes are mediated by subsurface microorganisms. To compare the bacterial and archaeal community compositions between groundwater groups, we performed both clone libraries construction and Illumina MiSeq sequencing of 16S rRNA genes from groundwater samples in springs or wells (Table S2). In a previous study, soil microbial community analysis using two different sequencing platform (i.e., clone libraries vs. pyrosequencing on a Roche 454 FLX machine) showed some similar results between methods, but also indicated critical differences in microbial community compositions (Jones et al., 2009).

Proteobacteria was the most abundant taxa in all water samples regardless of sequencing platform used, but bacterial communities followed by Proteobacteria were different between the methods (Table S3). In addition, although groundwater samples were extracted under the same conditions and analyzed using Illumina MiSeq platform, microbial community composition was different depending on primer sets used, which was consistent with the previous observation by Tremblay et al. (2015). Before MiSeq run, we performed RDP's probe match for each primer to evaluate the bias associated with each primer pair. Among bacteria, primer 341F (0.95 fraction) covered slightly more than primer 515F (0.93 fraction) (Table S4) and primer 341F generated more bacterial sequences than 515F in all samples (Table S2).

In general, analyses of MiSeq sequencing and clone library data yielded similar results. However, Illumina MiSeq sequencing revealed the greater diversity compared to clone library sequencing (Table S2). In addition, the MiSeq sequencing with 341F-805R yielded greater sequence coverage up to 10-fold depth than that with 515F-805R′. Therefore, we discussed the results of bacterial diversity and taxonomic assignments using the MiSeq sequencing with 341F-805R hereafter.

Analysis of rarefaction curves (plots of the number of species as a function of the number of individuals sampled) indicated the observed species at 97% similarity cut-off. The steep slope in all cases suggested that some fraction of the species diversity remains to be discovered (Figure 2). The rarefaction curves showed that the bacterial diversity was different between groundwater types. A rarefaction curve of groundwater samples influenced by high CO₂ concentrations revealed lower diversity, suggesting that elevated CO₂ level can increase the selective pressure for specific bacterial community (Figures 2A,B). Bacterial diversity indices (observed species) significantly increased from group II → group I → group III (Figure 2B). DPW-8, BG, and DPW-2 having high CO₂ indicated very low number of observed species (<1,100) (Table 2A). Chao1 and Shannon indices calculated based on OTUs also showed relatively low diversity in DPW-8 (620.4 and 5.0, respectively), BG (1160.6 and 7.5), and DPW-2 (790 and 6.5). Meanwhile, DPW-6 and DPW-7 in low-CO₂ control group III showed a very high number of observed species (>4,000). The Chao1 and Shannon indices of DPW-6 (4,499.6 and 9.5, respectively) and DPW7 (4,096.7 and 8.1) also indicated a relatively high diversity.

To confirm the variation of the species composition in this area, beta diversity calculated based on UniFrac analysis was compared between groundwater samples. PCoA plots based on beta diversity displayed three clearly distinguishable groups, which correspond to differences in water quality for bacteria (Figure 2E). In a complimentary analysis, NMDS analysis using top 100 bacterial genera (Table S5) results in the same clustering together (Figure 3A). There were significant differences (r² = 0.73, P < 0.01) between bacterial communities in groundwater according to their group as determined by PERMANOVA analysis based on Bray-Curtis distance matrix (Anderson, 2001). Environmental variables were overlaid to the plot to confirm whether geochemical factors affect community structure (Figure 3A). Among factors we tested, TDS, alkalinity, Sr, Li, Cs, and CH₄ showed some significant relationships with community compositions (0.54 < r² < 0.70, P < 0.1). The microbial community in group II was closely related to TDS, alkalinity, Sr, Li, Cs, and CH₄, while that in group III was correlated with relatively high pH. In addition, CO₂ showed close relationships with group I and group II suggesting that CO₂ may control bacterial community compositions in these groups.

Bacterial analysis by MiSeq sequencing with 341F-805R showed a total of 36 phyla as the following order of average relative abundances; Proteobacteria (40.2%) > Unclassified Bacteria (11.4%) > Bacteriodetes (9.4%) > Candidate division OD1 (7.7%) > Actinobacteria (6.4%) > Candidate division TM7 (5.4%) > Nitrospirae (3.3%) > Firmicutes (2.8%) > Cyanobacteria (2.6%) > Chloroflexi (1.8%) (Figure 4 and Table S3C). Proteobacteria was the most abundant phylum in all samples, but ranged broadly between 28.9 and 52.3%. Substantial variation of bacterial community compositions at the phylum level between water samples is not astonishing because different bacterial communities inhabit different environments. However, there are common taxa in specific groundwater group samples. For example, Firmicutes were dominant in group I (5.0 ± 4.3%) and Planctomycetes were abundant in group III (3.7 ± 0.9%) compared to other group samples, while Nitrospirae (9.2 ± 0.1%), Chloroflexi (4.0 ± 1.6%), and Chlorobi (3.0 ± 0.8%) were predominant in group II groundwater samples (Figures 3A, 4, Tables S3C, S6A–F). These results suggest that the specific bacterial phyla linked closely to the specific groundwater quality.

The predominance of the phylum Nitrospirae in group II was notable (Table S3C). Nitrospiraceae is the only established family in this phylum and contains three genera including Nitrospira (chemolithoautotrophic aerobic nitrite oxidizer), Leptospirillum (chemolithoautotrophic aerobic and acidophilic ferrous iron oxidizers), and Thermodesulfovibrio (anaerobic, thermophilic, chemoorganoheterotrophic or hydrogenotrophic sulfate reducers) (Table S6C). Although Nitrospiraceae has been isolated from and observed in ocean water, freshwater, aquarium water, marine sediment, soils, wastewater treatment plants, and a corroded iron pipe in a heating system (Watson et al., 1986; Ehrich et al., 1995; Daims et al., 2001; Altmann et al., 2003; Daims, 2014), their predominance in the group II groundwater samples (i.e., not thermophilic temperature, anaerobic, no nitrate) is enigmatic. The order Ignavibacteriales of phylum Chlorobi was also notable in group II (DPW-8; 3.5%, BG; 2.3%, respectively) (Table S6E). The family members in the order Ignavibacteriales are known to play important role in nitrogen cycle (Tian et al., 2015) including ammonification, nitrogen fixation and nitrification, and sulfur-oxidizing autotrophic denitrifying (Zhang et al., 2015).

Relatively high abundance of Chloroflexi, particularly Anaerolineaceae, was observed in group II (DPW-8; 2.2%, BG; 3.6%, respectively) (Table S6D). Anaerolineaceae has been reported as the most dominant taxa syntrophically associated with methanogenic archaea after long term incubation of methanogenic alkanes-degrader (Liang et al., 2015), and also belonged to dominant community in microbial electrolysis anaerobic digestion reactor system for accelerated methane production (Liu et al., 2016). The presence of Anaerolineaceae with methanogenic archaea in group II of our study are also consistent with those studies.

Group I groundwater samples (low pH, low TDS, and high CO₂) showed relatively high abundance of family Comamonadaceae (15, 8.7, and 9.6% in DPS-2, DPW-1, and DPW-2, respectively) (Table S6F). The result is consistent with a field-scale CO₂ geosequestration study that after CO₂ injection the relative abundance of Comamonadaceae increased with changes in the associated variables (e.g., pH, temperature and salinity) (Mu et al., 2014).

For archaeal community analysis, we also compared the sequencing results from three different methods (i.e., clone library with 109F-958R and MiSeq with 515F-805R′ or 341F-805R). Like bacterial analysis, we assessed the compatibility of two different primer sets for detecting archaeal sequences. The result of RDP's probe match showed that 515F covered 0.56 fraction of archaeal sequences but 341F did not cover at all (Table S4). However, 341F-805R indeed revealed a better depth of sequences than 515F-805R′ in this study (Table S2B). Table 2B summarizes the number of observed species and various diversity indices of archaea. Although the depth of sequences was reasonable for determining archaeal diversity, the sequencing result even with 341F-805R did not show the accurate information of archaeal community compositions because the most reads were unassigned archaeal OTUs in both reference databases, RDP and SILVA. The taxonomic assignment was better in clone library analysis than in sequencing either with 515F-805R′ or 341F-805R primer sets. Therefore, we discussed hereafter archaeal diversity based on the sequencing results with 341F-805R, but archaeal taxonomy using the clone library sequencing results with 109F-958R.

Similar to bacterial diversity, the alpha diversity slopes of archaea in the presence of high CO₂ were significantly lower than those in low CO₂ water. Archaeal diversity indices (observed species) also significantly increased from group II → group I → group III. The species diversity in partial pressure with no or relatively low concentrations of CO₂ was similar each other, but the presence of high CO₂ (>3,200 mg L⁻¹) decreased the number of observed species (Figures 2C,D).

The species richness of DPW-8 and BG in group II and DPW-2 in group I was low, which might be correlated with high CO₂ concentrations. DPW8 and BG showed low diversity indices of observed species, Chao1, and ACE (Table 2B). Unlike bacterial results, DPW-2 showed the lowest values of observed species (69), Chao1 (69.7), and ACE (72.2).

The beta diversity analysis for better understanding variation of species composition in archaeal community showed the three separated clusters in 3D-PCoA plot (Figure 2F). Similar to the bacterial community structure, DPW-8 and BG in group II and DPW-6 and DPW-7 in group III were very closely distanced, but DPS-2, DPW-1, and DPW-2 in group I were dispersed. This result suggests that microbial diversity in springs and wells is strongly affected by the level of CO₂ concentrations.

Archaeal community analysis by clone library construction with 109F-958R showed a total of 6 phyla. The result of NMDS from archaeal data (Figure 3B) revealed that relative abundance of top 16 archaeal genera (Table S7) across samples was linked distinctly to each groundwater group (r² = 0.53, P = 0.11). In addition, archaeal community compositions showed significant correlation with several environmental variables including alkalinity, Li, Rb, and Cs (r² > 0.89, P < 0.01) and TDS, SiO₂, NO3-, Sr, Cr, and CH₄ (r² > 0.82, P < 0.05).

Taxonomic distribution of archaeal community was distinct in each group (Figure 5 and Table S8). Euryarchaeota were dominated in group II groundwater samples, while Thaumarchaeota containing ammonia-oxidizing archaea (AOA) were predominant (73–88%) in group I and III groundwater samples. Ammonia-oxidizing bacteria (AOB) (Jiang and Bakken, 1999) and AOA (Könneke et al., 2005) are known to be able to oxidize NH₃ to NO2-. High abundance of AOA (e.g., Nitrososphaera and Nitrosopumilus) suggests that nitrification occurs to some extent in these environments. Nitrososphaera as acidophilic or thermophilic AOA usually inhabit terrestrial environments (i.e., soil), whereas Nitrosopumilus as mesophilic AOA usually presented in aquatic environments (i.e., marine water, freshwater and wastewater) (Spang et al., 2010; Pester et al., 2012). In addition, Nitrobacter, a genus of well-known nitrite oxidizing bacteria (NOB) (i.e., NO2- → NO3-), was observed in group I and III (i.e., DPS-2:0.1%, DPW-1:0.8%, DPW-6:0.1% and DPW-7:0.1%) (Teske et al., 1994). The presence of NOB as well as AOA with substantial amount of nitrate only in group I and III indicate the nitrification (ammonia and nitrite oxidation) in this environment (i.e., NH₃ → NO2- → NO3-).

Groundwater samples with high CO₂ and CH₄ showed clear differences in archaeal community compositions at the phylum level (Figure 5A). High abundance of methanogens in archaeal community in DPW-8 and BG coupled with presence of CH₄ suggests that methanogenesis is actively ongoing in group II groundwaters. Relative abundance of Euryarchaeota related to methanogenesis in DPW-8 and BG is 95.8 and 78.8%, respectively (Figure 5A), which corresponds well to the presence of CH₄ in these samples (35 and 22 mg L⁻¹, respectively) (Table 1). Interestingly, genus level compositions of archaeal community in between DPW-8 and BG were very different. Most archaeal sequences in samples from DPW-8 grouped in Methanobacterium (93.8%) whereas in samples for BG, archaeal sequences grouped in Methanosarcina (50%), Methanothrix (13.5%), and Methanomassiliicoccus (11.5%).

There are three known metabolic pathways for methanogenesis (Whitman et al., 2006); (1) hydrogenotrophic (e.g., 4H₂ + CO₂ → 2H₂O + CH₄), (2) acetotrophic (e.g., CH₃COOH → CO₂ + CH₄), and (3) methylotrophic (e.g., 4CH₃OH → 3CH₄ + CO₂ + 2H₂O). Methanobacterium species are strict anaerobes and hydrogenotrophic methanogens capable of using H₂ (Schirmack et al., 2014). Meanwhile, Methanosarcina and Methanosaeta (“Methanothrix”) are two genera among methanogens capable of using acetate to produce methane via acetoclastic methanogenesis (Zinder et al., 1984; Jetten et al., 1992; Kendall and Boone, 2006). Therefore, differences in the methanogenic archaeal community between DPW-8 and BG suggests that the methane in samples from each site may have been generated by different methanogenic pathways. Given that groundwater samples from each site had similar major ion chemistry and CO₂ content, we hypothesize that differences in methanogenic communities may be related to differences in DOC concentration in BG (4.7 mg L⁻¹) compared to that in DPW-8 (1.7 mg L⁻¹).

Our findings imply that microbial reactions in the shallow aquifer are influencing the carbon budget of the system in a variety of ways. Reaction of CO₂ with carbonate and silicate minerals, as discussed previously, generates carbonate alkalinity as well as divalent metals (e.g., Ca²⁺, Mg²⁺, Fe²⁺, and etc.) that can ultimately precipitate as carbonate minerals (e.g., calcite, dolomite, siderite, and etc.). Where those reactions occur, CO₂ would be trapped in mineral form, rather than emitted to the atmosphere. The microbial community may enhance mineral trapping through microbial generation of carbonate alkalinity. Many microbial reactions consume protons. As such, they drive CO₂ hydrolysis forward and thus convert CO₂ into bicarbonate. Microbial iron reduction, which may be generating the ferrous iron detected in samples from group II, is particularly effective at generating alkalinity because the reaction consumes several protons (Kirk et al., 2013). Microorganisms may also enhance carbon trapping by providing surfaces serve as nucleation sites for carbonate mineralization (Mitchell et al., 2010; Benzerara et al., 2011).

In addition to converting CO₂ into carbonate alkalinity or carbonate minerals, microorganisms can also transform the carbon in CO₂ into organic forms. Although the Wood-Ljungahl pathway isn't widespread within the Chloroflexi, two Chloroflexi species have been shown to use the Wood–Ljungdahl pathway for CO₂ fixation and play an important role in sediment carbon cycling (Hug et al., 2013). Therefore, the Chloroflexi may be helping to trap CO₂ as organic carbon in the subsurface at the site. In contrast, reduction of CO₂ to methane by hydrogenotrophic methanogens could result in loss of carbon as methane from the shallow subsurface. A relatively high abundance of Euryarchaeota associated with methanogenesis was observed in group II samples that also contained significant levels of CH₄. Methane is a potent greenhouse gas and its global warming potential is considered 25 times higher than that of CO₂ (Oppermann et al., 2010). Thus understanding microbial controls and developing engineered controls on its release from deep subsurface environments to the atmosphere is an important research need.

Alterations in microbial community structure have been observed during CO₂ storage in rhizosphere (Grayston et al., 1998) and soil (Rillig et al., 1997; Mayr et al., 1999; Horz et al., 2004) and thus can be used as a sign of CO₂ leakage. The previous studies of field-scale CO₂ geosequestration showed that the relative abundance of Comamonadaceae increased early after CO₂ injection with decreases in pH and mobilization in minor elements, such as Fe and Mn, temperature and salinity) (Mu et al., 2014; O'Mullan et al., 2015). Our study also showed clear relationships between microbial community compositions and CO₂-rich waters with relatively short-term or long-term water-rock interactions after the influx of CO₂ gases. As shown in group I case, the detection of Comamonadaceae in groundwater samples collected from near CO₂ storage sites would be an indication of the location influenced by CO₂ leakage with short-term water-rock reaction. The observation of methanogens, Nitrospirae, and Anaerolineaceae as shown in group II case may be an indication of the sites undergone both more CO₂ leakage and long-term water-rock reaction and biogeochemical reaction. Also, the decrease in bacterial and archaeal diversity in groundwater with increase in CO₂ would be an indication that CO₂ release from the storage sites is in progress and impacts groundwater microbial community. However, the differences in bacterial and archaeal diversity and compositions would be site specific. Therefore, long-term and periodical monitoring of microbial community in subsurface environments are necessary.