The Hetao Plain (40°10′–41°20′N, 106°10′–109°30′E) is one of the Mesozoic-Cenozoic faulted basins located in the western part of Inner Mongolia. It is along the northern bank of the Yellow River, and is bounded to the north by the Yin Mountains and to the west by the Wulanbuhe Desert, and with lacustrine plain in the central part (Figures 1A,B). Groundwater in the Hetao Plain is discharged mainly by vertical evapotranspiration and pumping and recharged by irrigation. Our case study was carried out in Hangjinhouqi County (Figures 1B,C) where endemic arseniasis is most serious in the Hetao Plain, with about 76 thousands local residents exposed under As-rich groundwater (Deng, 2008). In Hangjinhouqi County, four main drainage channels (three SE–NW and one SW–NE orientated) have been used to drain redundant irrigational water and to lower the local groundwater table (Figure 1C).

Groundwater samples along the agricultural drainage channels were collected from 17 domestic wells of seven arseniasis-affected villages in Hangjinhouqi County, of which sample IC12 and IC16 had been collected for the study of methanogens in high As groundwater in a previous study (Wang et al., 2015). Sampling locations were shown in Figure 1. Wells were pumped to get stabilized oxidation-reduction potential (ORP) values (commonly 10–15 min) prior to groundwater sampling. Microbial samples were collected by filtering of approximately 10 L fresh groundwater through 0.2-μm filters (Millipore), the filters were immediately packed into 50 mL sterile tubes, and then frozen in dry ice. All samples were transported to the laboratory on dry ice and then maintained at -80°C until further analysis. The geochemical parameters including temperature, pH, conductivity (COND), and ORP, were measured in situ using a multiple parameter water quality meter (Horiba, Japan). H₂S, NH4+, and Fe(II) were determined using a portable spectrophotometer (HACH, DR890) according to the manufacture’s protocols.

Sampling methods for ions and As speciation analyses were from Deng (2008). Briefly, groundwater samples were filtered through 0.45 μm mixed cellulose ester membranes and then the filtrates for cation analysis were acidified to pH < 2. Water samples were filtered through a disposable syringe joined with an anion exchange cartridge (Supelco, USA) to separate soluble arsenite [As(III)] and arsenate [As(V)] species, which were subsequently measured using liquid chromatography-atomic fluorescence spectrometry (LC-AFS-9700, Haiguang, China). Measurement of anions including NO3– and SO42– was made with ion chromatography (DX-120, Dionex, USA). Major cations were determined with ICP-AES (IRIS Intrepid II XSP). Iron species were determined by the 1,10-phenanthroline-based assay which was described in our previous study (Li et al., 2013). Total organic carbon (TOC) was measured using a TOC analyzer (Vario MICRO cube, Elemental, Germany). All samples were run in triplicate and then averaged.

For the analysis of microbial community composition, the filtered membranes were cut into small pieces under aseptic conditions and used for DNA extraction by the FastDNA SPIN Kit for Soil (MP Bio, USA) with a final elution in 80 μL deionized water. The universal forward primer 515 (5′-GTG CCA GCM GCC GCG GTA A-3′) and reverse primer 806 (5′-GGA CTA CCA GGG TAT CTA AT-3′) barcoded with a 12-base Golay code, which could cover a wide diversity of both archaea and bacteria (Bates et al., 2011), were used to amplify the microbial 16S rRNA gene V4 region. The 25 μL PCR reaction mix consisted of 2.5 μL of Takara 10× Ex Taq Buffer, 2 μL of dNTP Mix (2.5 mM each), 0.7 μL of Roche bovine serum albumin (20 mg/mL), 0.125 μL of Ex Taq DNA polymerase (Takara, 5 U/μL), 0.5 μL of 10 μM primer 515F, 0.5 μL of 10 μM barcode primer 806R, 1 μL of template DNA, and 17.675 μL of PCR grade water. The PCR amplification condition was: initial denaturation at 95°C for 3 min, followed by 25 cycles of 95°C for 30 s, 50°C for 30 s, and 72°C for 45 s, and a final extension at 72°C for 10 min. All samples were run in triplicate and then pooled together. Successful amplifications were confirmed by E-Gel (Invitrogen, USA) electrophoresis. PCR products were purified using the Agencourt AMPure XP PCR purification system (Beckman Coulter, Brea, CA, USA). The purified amplicons were quantified using the Qubit dsDNA HS assay and the size of the amplicons was determined using a Bioanalyzer with Agilent DNA 1000 chips (Agilent Technologies, Santa Clara, CA, USA). Amplicons were pooled (25 ng per sample) and the mixed suspensions were then subjected to paired-end (PE) sequencing by an Illumina MiSeq 2000 instrument at the Yale Center for Genome Analysis.

Principal component analysis (PCA) between As concentrations and other geochemical parameters were conducted using CANOCO for windows version 4.5. Hierarchical cluster analysis was performed with the “vegan” package in RStudio¹ using unweighted pair-group method with arithmetic means (UPGMA) and Euclidean distance measure. SPSS 20 software was used to analyze the significant differences (p ≤ 0.05) among different geochemical characteristics and microbial diversity.

The PE reads generated from Illumina MiSeq sequencing were assembled using FLASH (Fast Length Adjustment of Short reads) according to index sequence (Magoč and Salzberg, 2011). All sequences with ambiguous base calls were discarded. The PE sequences were also discarded if a mismatch was observed in the assembly. The output data was further analyzed using QIIME (quantitative insights into microbial ecology) v1.7.0. Operational taxonomic units (OTU) picking was performed using a “closed-reference” protocol. The most abundant sequence in each cluster was chosen to represent the cluster. All representative sequences were aligned with PyNAST. Chimeric sequences were identified with Chimera Slayer and excluded. Sequences were clustered against the Greengenes database² at 97% identity using uclust. Following assignment, 12,000 successfully assigned sequences from each sample were chosen at random to allow for downstream analyzes and cross-sample comparison. Alpha diversity analysis involving rarefaction curves, Chao1 and Shannon diversity (representing estimated OTU numbers and microbial diversities in the community, respectively) indices and beta diversity analysis were generated using OTU picking result in QIIME software.

Multivariate community analysis including principal coordinate analysis (PCoA) and co-inertia analysis was performed with the ade4 package (Chessel et al., 2004; Dray and Dufour, 2007) within R programming environment³ using normalized OTU data. Correlated geochemical variables were removed from the analyzed table. Besides, some geochemical variables were made Log₁₀ transformation to reach normal distribution before analysis. SIMPER (similarity percentage) analysis (Clarke and Warwick, 1994) was conducted to rank the top 10 OTUs that contributed to the dissimilarity index between high As groundwater samples and low As groundwater samples. The DNA sequences were deposited to the Short Read Archive database at NCBI (Accession number: SRR4996298).

Arsenic concentrations of groundwater in Hangjinhouqi County were generally elevated along the agricultural drainage channels in the front of Yin Mountains (Figure 1; Table 1). This was consistent with the results by He (2010) that As in soils or sediments could be leached into groundwater with the irrigation return flow. The 17 groundwater samples were divided into two groups with the distinctive geochemical characteristics according to the PCA (Figure 2A) and Hierarchical clustering results (Figure 2B). The first group (low As group) contained six groundwater samples which were characterized with low concentrations of As (As < 10 μg/L) and relatively low TOC and NH4+. The other group (high As group) included 11 samples with high concentrations of As (As > 10 μg/L), TOC and NH4+. Arsenic concentrations of those high As groundwater samples were from 50 to 1089 μg/L and As(III) was the dominant species in most of high As groundwaters, with ratios of As(III) to As_Tot ranging from 0.04 to 0.93 (Table 1). The concentrations of TOC lay in the range 1.9–18.4 mg/L in high As groundwater samples and 0–4.3 mg/L in low As samples. NH4+ concentrations varied from 1.94 to 6.68 mg/L in high As groundwater samples and from 0 to 3.87 mg/L in low As groundwater samples.

Arsenic concentrations of groundwater samples showed positive correlations with NH4+ and TOC (Table 2; Supplementary Figure S1). SO42– and Fe(II) concentrations in this study ranged from 35.41 to 1458.10 mg/L and 0.01 to 1.25 mg/L, respectively, showing no obvious correlations with As concentrations (Table 2). These results were different from those in our previous study of samples from strongly reducing areas in Hetao Plain that As concentrations showed negative correlations with the concentrations of SO42– and positive correlations with the concentrations of Fe(II) (Li et al., 2015). Considering the positive correlation among As, NH4+, and TOC, these results suggest that groundwater samples collected along the drainage channels might present distinct biogeochemical characteristics and redox properties. Agricultural activities such as N-fertilizer application and irrigation might also play important roles in As release and mobilization (Busbee et al., 2009; Mayorga et al., 2013). The reasonable explanation might be that the irrigation and drainage processes could leach As and dissolved organic matter from surficial sediments to the underlying aquifer, which caused the biogeochemical changes (He, 2010).

The microbial community diversities and phylogenetic structure of the 17 groundwater samples were analyzed by Illumina sequencing. A total of 509,554 full length V4 tags were obtained after removal of low quality and chimeric sequences and 24,938 OTUs were obtained at a distance of 0.03. A variety of taxa were observed at the 97% OTU level, with 329–2823 observed OTUs and 0.64–0.84 coverage values (Table 3). The abundance of OTUs was observed to be higher in the low As groundwater samples than in high As samples (Table 3). The rarefaction curves generated by QIIME were normalized to the minimum number (17,000) of sequences. The Chao1 and Shannon diversity ranged from 465.26 to 3420.42 and 4.08 to 9.32, respectively (Table 3), showing higher values in low As groundwater samples (Table 3). These results implied that high As groundwater samples presented lower microbial diversities which was consistent with our result of the 16S rRNA gene clone library (Li et al., 2013). Compared with groundwater samples from strongly reducing areas of Hetao Plain (Li et al., 2015), samples collected along the agricultural drainage channels showed much higher microbial diversities, implying that irrigation and drainage activities might affect local geochemical environment and thus enhance the diversity of microorganisms.

Microbial community compositions at phylum level showed no distinct difference between high and low As groundwater samples (Supplementary Figure S2). Illumina results of the 17 groundwater samples indicated that the bacterial community was composed of 62 phyla (percentages lower than 0.01% were not calculated). Of these, only nine phyla dominated each community (Supplementary Figure S2), among which Proteobacteria was the most dominant group (32.02–86.50%), followed by Firmicutes (0.16–18.48%), Actinobacteria (0.34–12.08%), and Nitrospirae (0–22.64%). Within phylum Proteobacteria, Gammaproteobacteria (6–83%), and Betaproteobacteria (5–39%) were the most abundant (Supplementary Figure S2). The Alphaproteobacteria and Deltaproteobacteria accounted for 0–17% and 0–23% of the total Proteobacteria, respectively. The class Epsilonproteobacteria displayed a low rate of occurrence and was only abundant in two high As groundwater samples IC02 and IC14 (Supplementary Figure S2). Archaeal populations accounted only for 0.01–8.69% of the microbial populations. All detected archaeal populations belonged to the phyla Euryarchaeota and Crenarchaeota which represented by the classes of Parvarchaea, Methanomicrobia, and Methanobacteria. Members of the Methanomicrobia and Methanobacteria retrieved at high As groundwater samples were mostly affiliated with methane-producing archaea Methanocorpusculum, Methanospirillum and SAGMEG-1 which have previously been detected in marine sediments and sludge (Iino et al., 2009; Barbier et al., 2012). These methanogenic populations might provide even stronger reducing conditions and thereby accelerate As release in groundwater aquifers in the Hetao Plain (Wang et al., 2015). However, irrigation and drainage activities in agricultural areas could bring more oxygen and the resulting oxidizing conditions might inhibit the activities of methanogenic populations.

At the genus level, the top 10 taxonomic OTUs that contributed to the dissimilarities between high and low As groundwater samples were identified by SIMPER analysis (Table 4). The average abundances of Alishewanella, Psychrobacter, Pseudomonadaceae, Methylotenera, Comamonadaceae, and Crenothrix were elevated in high As groundwater samples than in low As samples, while Rheinheimera and unidentified OP3 presented higher abundance in low As groundwater samples than high As samples (Figure 3B; Table 4). Although Acinetobacter and Pseudomonas had been reported to be correlated with As metabolism (Anderson and Cook, 2004; Chang et al., 2009; Freikowski et al., 2010), these two populations appeared with slightly higher percentages in low As groundwater samples but dominated in both high (0.2–28.2% and 0.1–13.6%, respectively) and low (0.2–41.4% and 0.2–37.3%, respectively) As samples in the present study. This might be due to the wide occurrence of these two kinds of bacteria in nature, and some of their isolates are capable of tolerating high concentrations of As. These results were mostly consistent with those of our previous studies in the strongly reducing area conducted with traditional sequencing and 454 pyrosequencing methods (Li et al., 2013, 2015). Predominant populations in high As groundwater samples such as Alishewanella and Psychrobacter were previously isolated from industrial effluent or groundwater aquifers, showing the capability of hyper arsenite tolerance or arsenate reduction (Liao et al., 2011; Jain et al., 2014). In some recent studies, species of Alishewanella were also identified as denitrifiers (Kolekar et al., 2013; Liu et al., 2015). Methylotenera species were reported as facultative methylotrophs and could be involved in the carbon and nitrogen metabolic pathways (Mustakhimov et al., 2013; Chistoserdova, 2014). The sample with the highest concentration of TOC presented the highest abundance of Methylotenera. Crenothrix was documented as one kind of iron oxidizing bacteria in groundwater (Kang and Liu, 2013; Thapa Chhetri et al., 2014). Besides, Crenothrix was also reported as a methane oxidizer with a unique methane monooxygenase (Stoecker et al., 2006). These dominant populations in high As groundwater indicated that As release in groundwater along agricultural drainage channels was correlated with microbial carbon, nitrogen, and iron reactions. Different from our previous studies with samples from the natural reducing area (Li et al., 2013, 2015), many oxidizing microbial populations such as Methylotenera and Crenothrix were found in the present study, indicating the possible influence of geochemical difference on microbial composition. The dominant populations in this study were also different from those in high As groundwater aquifer in Bangladesh (such as Hydrogenophaga and Acidovorax; Sutton et al., 2009), West Bengal (such as Agrobacterium–Rhizobium, Ochrobactrum, Anoxybacillus, and Paenibacillus; Sarkar et al., 2015), and Vietnam (Lawati et al., 2012). These community differences might be influenced by different geochemical conditions such as ORP (Sutton et al., 2009) and dissolved organic carbon (Lawati et al., 2012).

PCoA plots were used to visualize the relationships among microbial communities using default beta diversity metrics of unweighted Unifrac. Geochemical parameters including pH, ORP, As, Fe, SO42–, NH4+, and TOC were used to calculate the possible influence on microbial differences. Results showed that microbial communities could be generally divided into two groups according to different As concentrations (Figure 4). A similar result was also obtained with UPGMA cluster tree analysis which was based on Bray-Curtis dissimilarity at the 97% similarity OTU level (Figure 3A). These results implied that As was an important factor shaping the microbial community structure in groundwater of agricultural irrigation area.

Co-inertia analysis was carried out in order to reveal the relative importance of geochemical vectors that affect microbial community structures. The arrow head and end represented the position of sample according to the geochemical ordination and microbial ordination, respectively (Figure 5). As shown in Figure 5B, the first two axes explained 81% of the total environmental variability. Results revealed that high As groundwater and low As groundwater samples were significantly different (p value < 0.001) in community structures (Figure 5A) and geochemical variables had a significant influence on their microbial community compositions. In addition to As concentration, geochemical variables such as TOC, NH4+, Fe, and ORP might also affect the composition (Figure 5C). These important geochemical factors shaping the microbial community structure were different from those of strong reducing area in the Hetao Plain which were characterized with high concentrations of Fe(II), H₂S, and SO42–.