Preferential flow paths shape the structure of bacterial communities in a clayey till depth profile.

Preferential flow paths in subsurface soils serve as transport routes for water, dissolved organic matter and oxygen. Little is known about bacterial communities in flow paths or in subsoils below ∼4 m. We compared communities from preferential flow paths (biopores, fractures and sand lenses) with those in adjacent matrix sediments of clayey till from the plough layer to a depth of 6 m. 16S rRNA gene-targeted community analysis showed bacterial communities of greater abundance and diversity in flow paths than in matrix sediments at all depths. Deep fracture communities contained a higher relative abundance of aerobes and plant material decomposers like Nitrospirae, Acidobacteria and Planctomycetes than adjacent matrix sediments. Similarly, analyses of the relative abundances of archaeal amoA, nirK and dsrB genes indicated transition from aerobic to anaerobic nitrogen and sulphur cycling at greater depth in preferential flow paths than in matrix sediments. Preferential flow paths in the top 260 cm contained more indicator operational taxonomic units from the plough layer community than the matrix sediments. This study indicates that the availability of oxygen and organic matter and downward transport of bacteria shape bacterial communities in preferential flow paths, and suggests that their lifestyles differ from those of bacteria in matrix communities. © FEMS 2019.

INTRODUCTION

Soils display considerable heterogeneity and can be considered as several specific environments that harbour a large diversity of microorganisms (Fierer 2017). Whereas the factors shaping soil microbiome structure and function have been widely studied in topsoil, knowledge about bacterial communities in deeper soil is very limited (Sagova-Mareckova et al. 2016; Pereira et al. 2017). The few studies of subsoil communities have consistently reported a decrease in bacterial abundance and diversity in depths down to ∼1.8 m below ground surface (mbgs) (Agnelli et al. 2004; Will et al. 2010; Eilers et al. 2012). A decrease in diversity may imply an ecosystem relying on fewer functions and providing fewer services such as organic matter decomposition (Delgado-Baquerizo et al. 2016; Maron et al. 2018) or nutrient cycling (Philippot et al. 2013). However, the composition and functional potential of communities from layers below ∼1.8 mbgs (hereafter referred to as sediments) are still almost completely unknown.

In clayey tills, the presence of hydraulic preferential flow paths further increases heterogeneity and hence has an impact on bacterial diversity (Curd et al. 2018). Preferential flow paths may be biopores such as earthworm burrows and root channels (Jarvis 2007), tectonic fractures created by chemical and physical processes and, in some tills, sand lenses interspersed between the fractures. Biopores dominate down to 2 mbgs, while tectonic fractures become more prevalent at greater depths. Tectonic fractures exist down to 8 mbgs in various tills where they may be relatively abundant, with distances between fractures ranging from 0.1 to 2.5 m (Jørgensen, Klint and Kistrup 2003).

Biopores, fractures and sand lenses have greater hydraulic conductivity than the adjacent matrix sediment (Kessler et al. 2013), and in some cases may account for >94% of water transport (Jørgensen et al. 2002). Biopores, fractures and sand lenses may be interconnected, thereby serving as preferential pathways for the transport of water containing dissolved oxygen, nutrients, organic matter and even microorganisms (Rosenbom et al. 2008; Nielsen et al. 2010; Arnaud et al. 2015). However, the extent of this anticipated transport of microorganisms from topsoil to lower fractures remains unknown.

Plant exudates and material deposited by earthworms make biopores nutrient-rich hotspots for microbial activity (Kuzyakov and Blagodatskaya 2015), including carbon (C), nitrogen (N) and phosphorous (P) cycling (Uksa et al. 2014). Furthermore, they give rise to bacterial communities that are different from those in the matrix soil (Stromberger, Keith and Schmidt 2012; Uksa et al. 2015). To date, studies of bacterial communities in flow paths have been limited to biopores in the surface soil (0–1 m) (Uksa et al. 2015), possibly due to difficulties in sampling from deeper fractures. The connectivity of biopores and fractures suggests that even deep fractures may be hotspots for microbial activity compared with matrix sediments where there is much less water flow (Jørgensen et al. 2002). Hence, further analyses of both the composition and functional potential of bacterial communities from deeper fractures and sand lenses will be instrumental in identifying the factors which shape communities in these deep soil compartments, and in deciphering their role in biogeochemical cycling and their contribution to subsoil ecosystem processes (Seuradge, Oelbermann and Neufeld 2017).

The current study aimed at determining the composition and functional potential of bacterial communities in deep sediments and preferential flow paths in clayey tills, including biopores, fractures and sand lenses. Further, we investigated if bacterial communities in preferential flow paths were influenced more by bacteria in the plough layer than bacterial communities in matrix soils. The bacterial communities were characterised by 16S rRNA gene amplicon sequencing and their functional potential determined using a novel high-throughput qPCR method that targeted genes involved in the cycling of C, N, P and sulphur (S) (Zheng et al. 2018). It was hypothesised that preferential flow path communities have a different composition and greater diversity, both taxonomically and functionally, than matrix communities, due to their exposure to organic matter, nutrients and oxygen transported along preferential flow paths.

MATERIAL AND METHODS

Site and sampling

Sampling was carried out at an agricultural field in Lund, Denmark (55°14’49’’N, 12°17’24’’E) over a three-day period in September 2016. The Lund field is a clayey till which was deposited on top of fractured chalk during two advances during the Late Weichselian glaciation. The site is a part of the Danish Pesticide Leaching Assessment Programme. The lithology below the plough layer is characterised as an oxidised sandy till rich in biopores, mainly in the form of earthworm burrows penetrating to depths of 1.5 mbgs (Fig. 1). Vertical fractures were also present, being more dominant at depths below 1 m. The oxidised till became more clay-rich, with chalk clasts at depths below 2.5 m and fractures penetrating the reduced clay below. Small, dispersed sand lenses were found in both the oxidised and reduced till. The matrix sediments were dominated by silt material (40–50%) and fine to medium sand (30–50%). Clay increased from 5% in the plough layer to 15% at 4–6 mbgs. Total organic carbon was 1% in the plough layer and decreased to <0.2% at 0.5 mbgs and below (Harder et al. in press).

Figure 1.

Conceptual model of the excavation with the different sediment types and pictures from the sampling. Only major fractures are plotted. Arrows indicate the zones from where samples were taken. Modified from Harder et al. (in press).

A 6-m deep multiple bench excavation permitted sampling of soil and sediments at different depths. Sampling began immediately after the excavation was established. The walls of the excavation were covered by a plastic membrane between sampling days to avoid dry-out and contamination from rainfall. Samples were obtained from the four sides of the excavation in the plough layer at 20–30 cm below ground surface (cmbgs), earthworm burrows (biopores) and matrix soil at 70–120 cmbgs, fractures and matrix sediments at 200–260 cmbgs, 300–350 cmbgs and 450–480 cmbgs, and sand lens and matrix sediments at 550–600 cmbgs (Fig. 1, Table 1). The fractures and biopores were identified by visual inspection due to differences in colour between these structures and the corresponding matrix sediments (Fig. 1). No biopores were visible in the plough layer due to agricultural tillage.

Table 1.

Sampling depths and sample types. Sampling depth is shown as the mean (sampling interval). Preferential flow paths are abbreviated below as ‘pref. flow paths’.

Sample name	Depth [cmbgs]	Compartment	Lithology
PL25	25 (20–30)	Plough layer	Organic rich clay, no biopores
M95	95 (70–120)	Matrix	Reddish sandy till
F95	95 (70–120)	Pref. flow paths	Biopores in reddish sandy till
M230	230 (200–260)	Matrix	Reddish clay till
F230	230 (200–260)	Pref. flow paths	Grey fractures in reddish clay till
M325	325 (300–350)	Matrix	Reddish clay till, with chalk clasts
F325	325 (300–350)	Pref. flow paths	Reddish fractures in reddish clay till
M465	465 (450–480)	Matrix	Reduced greyish clay till, with chalk clasts
F465	465 (450–480)	Pref. flow paths	Reddish fractures penetrating reduced clay till
M575	575 (550–600)	Matrix	Reduced greyish clay till, with chalk clasts
SL575	575 (550–600)	Pref. flow paths	Reddish sand lens in reduced clay till

To avoid contamination, the outermost layer of the walls of the excavation was scraped off (min. 5 cm) before sampling. It was possible to obtain small amounts of sample material from the biopores, fracture linings, sand lens and adjacent matrix sediments using a sterile spoon. For each location, 10–15 samples were combined and passed through a 2-mm sieve to remove stones and roots. Subsequently, three subsamples (250 mg) were taken randomly from each combined sample and stored at −20°C until DNA extraction.

DNA extraction and sequencing of the bacterial 16S rRNA gene

DNA was extracted from 0.25 g samples with the Mobio PowerLyzer PowerSoil DNA Isolation Kit (Mobio, Carlsbad, California, USA) using the manufacturer's protocol. Salmon sperm G2 was added to the samples prior to extraction to increase extraction efficiency (Paulin et al. 2013). The DNA concentrations, including G2, obtained were in the range of 4.4–12 ng (g wet wt soil)−1, as determined using a NanoDrop UV-VIS Spectrophotometer (Saveen & Werner AB, Limhamn, Sweden). Amplicon libraries were prepared using the barcoded primers 341F (5΄-CCTACGGGNGGCWGCAG-3΄) and 806R (5΄-GACTACHVGGGTATCTAATCC-3΄), targeting the hypervariable V3-V4 regions of bacterial 16S rRNA gene (Yu et al. 2005). The resulting libraries were submitted to an Illumina MiSeq (Macrogen, Seoul, South Korea) for sequencing with a read length of 2 × 300 bp.

Quantitative PCR of 16S rRNA gene

The 16S rRNA gene copy numbers in the samples were determined by quantitative PCR using the primers 341-F (5΄CCTACGGGAGGCAGCAG-3΄) and 518R (5΄-ATTACCGCGGCTGCTGG-΄3) (Muyzer, De Waal and Uitierlinden 1993). As a standard, E. coli K12 was grown in Luria-Bertani medium for one day at room temperature and added to the soil before DNA extraction. The cells were counted before addition to the soil. A 10-fold dilution series of the DNA extracted from the soil with added K12 cells was prepared, ranging from 5.9 × 101 to 5.9 × 106 16S rRNA gene copies per µl, and used as a standard. Each PCR reaction (20 µl) was performed using 2 µl (10 mg/ml) bovine serum albumin (Bioron, Ludwigshafen, Germany), 0.8 µl of both forward and reverse primer (10 pmol), 5.4 µl sterilised deionised water, 10 µl SYBRgreen Lo-ROX (PCR Biosystems, London, UK) and 1 µl DNA. The reaction conditions were as follows: 2 min at 95°C for enzyme activation, followed by 34 cycles of 10 s at 95°C, 30 s at 50°C and 1 min at 72°C. The final elongation step was at 72°C for 6 min, followed by a melting curve analysis starting at 55°C, and reaching 95°C by incremental increases of 0.5°C s−1. The qPCR was performed on a Biorad CFX Connect Real-Time System (Bio-Rad Laboratories, Inc., Hercules, California, USA).

Sequence analysis

The raw sequences had the barcodes trimmed off and the sequencing data quality was checked using FastQC (Andrews 2010). Reads were merged and primers trimmed off with Pandaseq (quality threshold value of 0.9) (Masella et al. 2012). The resulting sequences were 375–475 bp. A total of 2 977 212 sequences were obtained from the merging step. Further processing of the data was performed using QIIME version 1.9.1 (Caporaso et al. 2010). Operational taxonomic units (OTUs) were picked using only non-unique sequences with usearch from the Uparse pipeline (Edgar 2013) with a 97% identity level. Chimeras (0.9% of the sequences) were removed. The sequences were aligned against the Greengenes database (DeSantis et al. 2006) and sequences not matching the database (32.6%) were left out of the final OTU table. In total, 1 998 362 sequences were present in the 33 samples, varying between 50 722 and 74 477 sequences per sample. A total of 12 719 unique OTUs (97% identity level) were present in the samples. The sequences represented 10 295 bacterial OTUs from 64 phyla at the 97% identity level; 113 OTUs were classified as Archaea and 2311 OTUs were classified as unassigned. Singletons were removed before further processing to reduce noise (Edgar 2013). The raw sequences have been deposited as a NCBI Bioproject under accession number PRJNA453676.

High-throughput qPCR

Seventy-one functional genes involved in cycling of C, N, P and S as well as the 16S rRNA gene were quantified in parallel by a high-throughput qPCR (HT-qPCR), quantitative microbial element cycling (QMEC) method, developed by Zheng et al. (2018). We have included the primer sequences in Table S1. We refer to the article for more information about the tool, primer specificity, standards, etc. Briefly, DNA extracts were diluted to 50 ng µl−1 with sterilised water before setting up the HT-qPCR. The primers and samples were automatically dispensed into a chip and quantification was performed on a SmartChip real-time PCR system (WaferGen Biosystems, Fremont, California, USA) according to the manufacturer's instructions. The amplification protocol was as follows: 10 min at 95°C, followed by 40 cycles of 30 s at 95°C, 30 s at 58°C and 30 s at 72°C. The WaferGen software automatically generated melting curves. To ensure amplification accuracy, results with multiple melting peaks, <80% or >120% amplification efficiency or a threshold cycle (CT) > 31 were excluded. Relative gene copy numbers were calculated using the equation: gene copy number = (31-CT)/(10/3). Functional genes were normalised to the amount of 16S rRNA genes to obtain relative gene abundance (gene copy numbers of functional genes per 16S rRNA gene). Water samples were included as negative controls.

Data analyses

All data analyses were conducted in R v. 3.3.1 (The R foundation for Statistical Computing, Vienna, Austria). Bacterial richness and diversity indices (Chao1 and Shannon) were calculated using the phyloseq package (McMurdie and Holmes 2013) without rarefaction (McMurdie and Holmes 2014). Rarefaction curves were generated using the vegan package (Oksanen et al. 2018). Archaea was not included in any of the analyses. The abundance of 16S rRNA gene copy numbers and bacterial diversity indices in matrix sediments and preferential flow paths were compared with one-tailed paired t-tests. Samples from the plough layer were omitted because they did not belong to either of the above compartments. The differences between sediment types were not normally distributed for the Shannon index and no paired t-test was applied. 16S rRNA copy numbers were log10 transformed before the paired t-test. P values <0.05 were considered statistically significant. Bar plots were made using the phyloseq (McMurdie and Holmes 2013) and ggplot2 (Wickham 2009) packages. OTUs were grouped at phylum level, and phyla with a relative abundance <2% were merged into ‘Others’. Unassigned OTUs were grouped separately and visualised if the relative abundance was >2% of the community in a sample. Non-metric multidimensional scaling (NMDS) was plotted using Bray-Curtis dissimilarities for the bacterial OTUs and functional genes with the ‘ordinate’ function in the phyloseq package. For the analysis of the functional genes, a minimum of two out of three replicates had to be above the detection limit in at least three sediment types to be included in the analysis. Correlation between the OTUs and functional genes was tested using the Mantel test in vegan (Oksanen et al. 2018). Heatmaps were generated using the ampvis package (Albertsen et al. 2015). Indicator OTUs were defined as bacterial OTUs having a higher abundance in one sediment type than in any of the other sediment types (Dufrêne and Legendre 1997). OTUs with a minimum abundance of three in at least three sediment types were admitted into the analysis for indicator OTUs. Indicator OTUs were calculated using the ‘indval’ function from package labdsv (Roberts 2016) at a significance level of p ≤ 0.05 and IndVal > 0.4, using the three samples from each sediment type as replicates. The indicator value was set to >0.4 to include more than just the unique indicator OTUs (indicator value = 1) for each sediment type. Differences in relative abundance of indicator OTUs between sample types were tested using Welch's t-test.

The nestedness temperature T was calculated using the algorithm developed by Atmar and Patterson (1993). The temperature is quantified on a scale of 0–100, with a perfectly nested matrix having a temperature of 0°C. In this study, the nested rank order of samples was determined using BINMATNEST (Rodriguez-Girones and Santamaria 2006) with the following parameters: Null matrices= 100, PopSize= 30, TourSize= 7 and nbGen= 2000. This analysis establishes if there is significant trend of nestedness related to sediment type.

RESULTS

Bacterial abundance and diversity

Bacterial abundance in the plough layer was 4.93 ± 2.85 × 108 DNA copies (g wet wt soil)−1, increasing to 1.28 ± 0.98 × 109 DNA copies (g wet wt soil)−1 in the biopores (Fig. 2A). In the underlying fractures the bacterial abundance decreased ∼1000-fold with the lowest copy number, 9.86 ± 1.26 × 105 copies (g wet wt soil)−1, measured at 465 cmbgs. An 80-fold increase in bacterial abundance was observed from the deepest fractures to the sand lens at 575 cmbgs (7.63 ± 2.55 × 107 copies (g wet wt soil)−1). The variation with depth for bacterial abundances in the matrix sediments was comparable to the variation observed for preferential flow paths; however, at all depths the preferential flow paths showed significantly higher bacterial abundance than the corresponding matrix sediments based on one-tailed paired t-test (t = 5.04, df = 4, P = 0.0036).

Figure 2.

(A) Vertical profiles of 16S rRNA gene copy numbers (log gene copies (g wet wt soil)−1 ), (B) Chao1 richness and (C) Shannon diversity indices. Plots show the measured variables for sampled compartments and depths. The mean depths of the sampling sites are used to plot points against the y-axis. Points are mean values (n = 3) with standard error of the mean. The plots share the y-axis.

The diversity indices and the rarefaction curves showed that bacterial communities in the plough layer and biopores had the highest diversity (Fig. 2, Fig. S1). The diversity decreased with depth in the underlying fractures and sand lens. Significantly lower Chao1 indices were measured in the matrix sediments compared with the preferential flow path samples, although similar values were obtained at a depth of 325 cmbgs (one-tailed paired t-test, t = −3.28, df = 4, P = 0.015) . The bacterial diversity (Shannon index) showed a similar trend to the Chao1 richness, as higher indices were obtained in the preferential flow paths than in matrix sediments (Fig. 2C).

Variation in bacterial community structure

Differences in community composition in the samples were measured using NMDS analysis based on the Bray-Curtis dissimilarity index (Fig. 3). The ordination of 16S rRNA gene data showed that the community in the plough layer was closely related to the community in the biopores at 95 cmbgs. Furthermore, the community from the sand lens at 575 cmbgs clustered with the community in the matrix sediment at 95 cmbgs. Besides these relationships, the NMDS plot showed limited clustering of samples.

Figure 3.

Non-metric multidimensional scaling (NMDS) ordination of Bray-Curtis dissimilarity on 16S rRNA gene abundance (stress = 0.083).

Proteobacteria (mean relative abundance: 22%), Actinobacteria (14%), Bacteroidetes (13%), Firmicutes (12%) and Chloroflexi (10%) were the most abundant phyla across samples (Fig. 4, Fig. S2). However, the abundance of individual phyla varied considerably between samples, as also seen for bacterial groups at lower taxonomic levels (Figs. S3-S5). In the plough layer, as well as in the biopores and matrix sediments at 95 cmbgs, the most abundant phyla were Proteobacteria, Actinobacteria, Bacteroidetes, Acidobacteria and Firmicutes. The relative abundance of Proteobacteria peaked at 230 cmbgs in the fractures (25%), and at 575 cmbgs in the matrix sediments (31%). The abundance of Proteobacteria showed distinct differences between the matrix and fractures at depths below 260 cmbgs. At 325 cmbgs their abundance was higher in the fractures than in the adjacent matrix, while their relative abundance was higher in the matrix sediments at 465 and 575 cmbgs.

Figure 4.

Relative abundance of bacterial phyla based on sequencing data of the 16S rRNA gene. Phyla <2% are grouped in ‘Others’. OTUs not classified at phylum level are grouped in ‘Unassigned’.

In the plough layer and at 95 cmbgs (matrix and biopores), Alphaproteobacteria was the most dominant proteobacterial class while Beta-, Delta- and Gammaproteobacteria had comparable relative abundances of 4–6% (Fig. S3). However, with increasing depth, the differences between their abundances became notable. For example, in the fractures and matrix sediment at 465 cmbgs, Deltaproteobacteria and Betaproteobacteria were respectively the most abundant, while at 575 cmbgs, Gammaproteobacteria and Alphaproteobacteria dominated in the matrix sediment and sand lens, respectively. Betaproteobacteria also dominated in both fractures and matrix sediment at 230 cmbgs, with abundances of 12% and 17%, respectively. In the matrix sediment at 465 cmbgs, the betaproteobacterial iron oxidisers Gallionellales made up 12% of the entire community, but <1% of the community in the adjacent fractures (Fig. S4).

The Actinobacteria decreased in abundance with depth until 325 cmbgs, where abundances of 9% and 7% were observed in fractures and matrix sediments, respectively. In deeper layers, the abundance of Actinobacteria increased again, reaching a maximum in the matrix sediment at 575 cmbgs (41%), while in the adjacent sand lens this bacterial group only comprised 15% of the community (Fig. S2). Whereas the two main Actinobacteria classes, Actinobacteria and the strictly aerobic Thermoleophilia, had similar abundances at 230 cmbgs in the sediment profile, Thermoleophilia was the most abundant in the fractures and Actinobacteria in the matrix sediments at greater depths (Fig. S3).

The distributions of the Bacteroidetes and Firmicutes were substantially different from those of Actinobacteria and Proteobacteria. Regarding variation in relative abundance with depth, no general trend was observed for these phyla. The maximum abundance of the Firmicutes was observed in the matrix sediments at 325 cmbgs (34%), while they were almost absent in the matrix sediment and fractures at 230 cmbgs (Fig. S2). Among the Firmicutes, the facultative anaerobic Bacilli was the most dominant class in the plough layer, biopores, fractures and matrix sediment at 230 cmbgs, while obligate anaerobic Clostridia was the major class in communities in the matrix sediment at 95 cmbgs and in the deeper sediments (Fig. S3).

Bacteroidetes showed maximal relative abundance in the matrix sediment at 325 cmbgs, while both the matrix and the fracture communities at 230 cmbgs only consisted of 2–3% Bacteroidetes (Fig. S2). In the plough layer, biopores and fractures at 230 cmbgs, the Cytophagaceae and Chitinophagaceae families were relatively more abundant compared with the matrix sediments in the top 230 cm (Fig. S5). In contrast, Bacteroidaceae and Prevotellaceae were notably more abundant in the matrix sediment at 95 cmbgs. The fall in the relative abundance of Bacteroidetes at 230 cmbgs was followed by an increase and clear shift in dominating families for the deeper sediments, where the anaerobic Bacteroidaceae became the most abundant family.

The abundance of Acidobacteria was stable to a depth of 230 and 325 cmbgs for the matrix sediments and fractures, respectively. Below these depths, the Acidobacteria decreased and they constituted <2% of the communities in the matrix sediments at 465 and 575 cmbgs (Fig. S2). Among the Acidobacteria, subdivision 6 was the most dominant at most depths (Fig. S3).

The abundance of the chemolithoautrophic aerobic nitrite-oxidising Nitrospirae varied with depth. Its relative abundance showed low values in the plough layer, increased at 95 cmbgs, and peaked in the fractures and matrix sediment at 325 and 230 cmbgs, respectively. Below this depth, Nitrospirae almost disappeared from the matrix sediments, while high abundances of the phylum were still found in the fractures at 465 cmbgs. The variation in relative abundance of the Gemmatimonadetes, and to a certain extent the plant material decomposing Planctomycetes, resembled the variation recorded for Nitrospirae, although they were both found in lower abundances. Chloroflexi increased with depth to 230 cmbgs. At 95 cmbgs they were more abundant in the matrix sediments, but below this depth the opposite was observed. Verrucomicrobia was primarily detected in the upper sediments (plough layer and 95 cmbgs), but also in the sand lens at 575 cmbgs.

Relative abundance of genes involved in N and S cycling

The 71 functional genes' abundances involved in C, N, P and S cycling were obtained by HT-qPCR. NMDS analysis showed that the plough layer, biopores and matrix (95 cmbgs) and sand lens clustered, while the fractures and matrix sediments at 230 and 465 cmbgs were separated (Fig. 5). The dissimilarity between clustering of functional and 16S rRNA genes was verified by the Mantel test (R = 0.389, P < 0.001). For most C, P and S cycling genes, the relative abundance increased with depth (Figs. S6-S11), e.g. the manA gene, encoding mannose-6-phosphate isomerase, which increased from 0.0027 genes per 16S rRNA in the plough layer to 0.009 in the fractures at 465 cmbgs and the sand lens (Fig. S7).

Figure 5.

Non-metric multidimensional scaling (NMDS) ordination of Bray-Curtis dissimilarity on functional gene abundance from the HT-qPCR chip, where genes detected with low frequencies were removed (stress = 0.052).

The present study focused on the relative abundances of genes involved in N and S cycling. The amoA1 gene, encoding archaeal ammonia monooxygenase subunit A, had high relative abundances in the matrix sediment communities down to 230 cmbgs. For the preferential flow paths, they were detected to depths of 325 cmbgs and, interestingly, in the sand lens at 575 cmbgs (Fig. 6A). The relative abundance of amoA1 reached a maximum in the matrix sediment at 230 cmbgs and in the fractures at 325 cmbgs, respectively (Fig. S9). Other genes involved in different steps of ammonia oxidation and nitrification (bacterial amoA2, amoB and nxrA), were primarily detected in the top 95 cm of the sediment profile and in the sand lens (Fig. S9).

Figure 6.

(A) Relative abundance of archaeal amoA1, nirK1and nosZ1 genes per 16S rRNA gene and (B) relative abundance of bacterial soxYanddsrB genes per 16S rRNA gene. Error bars are standard error of the mean (n = 3).

For genes involved in reductive N metabolism, the nirK1 gene (encoding a nitrite reductase) had a lower relative abundance than the amoA1 gene down to 230 cmbgs for both preferential flow paths and matrix sediments (Fig. 6A). In contrast, the abundance of nirK1 was higher than amoA1 in the deeper sediments. The highest abundance of the nirK1 gene was found in the matrix sediment at 465 cmbgs, where other genes involved in reductive N processes, including nirK2, nirS1and nirS2 genes, also had optima (Fig. S9). However, the nosZ2 gene, encoding nitrous oxide reductase, and the nosZ1 and hzsB genes involved in denitrification and anammox processes, respectively, were not detected in the matrix sediment at 465 cmbgs (Fig. 6A, Fig. S9).

The relative abundance of the dsrB gene encoding the beta subunit of sulphite reductase (involved in the reduction of sulphite to sulphide) increased substantially at 325 cmbgs and reached a maximum in the matrix sediment at 465 cmbgs (Fig. 6B, Fig. S10). For comparison, this gene was less abundant in the fractures at this depth. Similarly, the abundance of the soxY gene involved in the S-oxidising pathway increased with depth, reaching a maximum in the matrix sediment at 325 cmbgs and in the deeper preferential flow paths at 465 cmbgs (fractures) and 575 cmbgs (sand lens) (Fig. 6B).

Indicator OTU analysis

The plough layer community had the highest number of indicator OTUs (445), while the sand lens and matrix sediments at 575 cmbgs had the lowest (48 and 26, respectively). No unique indicator OTUs, which were only present in one compartment, were detected in any of the compartments. While the biopore community consisted of 11% indicator OTUs from the plough layer, the matrix community at 95 cmbgs had <2% (Fig. 7), which was significantly lower (t-test, P = 0.006). Additionally, the preferential flow paths at 230 and 575 cmbgs contained significantly more indicator OTUs from the plough layer than the corresponding matrix sediments (t-test, P = 0.028 and P = 0.036, respectively). At 325 cmbgs, indicator OTUs from the plough layer did not show any differences in relative abundance between the matrix sediment and fractures.

Figure 7.

Mean abundance (%) of plough layer indicator OTUs in the different sediment types. * significant at P < 0.05; ** significant at P < 0.01.

The observation that indicator OTUs from the plough layer were more dominant in some preferential flow paths suggests a colonisation of the deeper communities by bacteria from the plough layer, possibly due to transport via preferential flow paths. To investigate whether there was a downward hierarchical organisation of the communities, we performed a nestedness analysis in order to determine whether the deep matrix and preferential flow path communities were subsets of the upper matrix and preferential flow path communities, respectively. The nestedness temperature was 28.3 (P < 0.001; null model 3). The test showed that the plough layer and biopore communities had the highest orders, and that the preferential flow path communities underneath were subsets of these communities (Fig. S12). The deep matrix communities had the lowest orders, indicating that they were subsets of the upper matrix sediment communities. Both the distribution of the indicator OTUs and the results of the nestedness analysis showed that plough layer communities influenced the biopore, fracture and sand lens communities more than they influenced the matrix sediment communities.

DISCUSSION

Bacterial abundance, richness and diversity decreased with depth both in preferential flow paths and adjacent matrix sediments, as observed previously for matrix sediments (Eilers et al. 2012; Sagova-Mareckova et al. 2016; Seuradge, Oelbermann and Neufeld 2017). Furthermore, bacterial abundance and richness were significantly higher in the biopores, fractures and sand lens compared with the adjacent matrix sediments. It is likely that the input of dissolved oxygen, nutrients and organic C decreases with depth, while the biopores and fractures, serving as preferential flow paths for water, provide these resources for bacterial growth from the surface (Rosenbom et al. 2008; Nielsen et al. 2010; Arnaud et al. 2015; Seuradge, Oelbermann and Neufeld 2017). Due to the general low porosity of clayey till, the transport of nutrients and organic C into the matrix may primarily be governed by diffusive processes, which are most likely too slow to support diverse microbial populations.

Soil organic C, as well as the availability of oxygen, may play an important role in the bacterial community composition in the current deep sediment layers (Lüdemann, Arth and Liesack 2000; Hansel et al. 2008). The distribution of several taxa seemed to be influenced by oxygen conditions. Hence, the aerobic chemolithoautotrophic nitrite oxidisers Nitrospirae were abundant in all samples from preferential flow paths, while they almost disappeared in matrix sediments at 325 cmbgs and below, most likely due to a transition from aerobic to anaerobic conditions. Further, these bacteria do not benefit from soil organic carbon in the top soil layers, and Nitrospirae has been shown to be suppressed by root-associated microbiota, which may explain its low relative abundance in the plough layer and biopores (Will et al. 2010).

In addition to the Nitrospirae, Acidobacteria subdivision 6 harbouring taxa, which can respire oxygen at low sub-atmospheric concentrations (Vieira et al. 2017; Eichorst et al. 2018) were more abundant in the preferential flow paths than in the matrix sediments at depths >230 cmbgs. Furthermore, they harbour the ability to metabolise complex organic molecules, and genomic analyses have suggested that some strains could also fix C (Eichorst et al. 2018). Hence, Acidobacteria subdivision 6 have the genetic potential to thrive in low oxygen conditions, which may occur in fractures, while having the ability to fix C if other C sources are depleted.

The distribution of lower taxa belonging to Bacteroidetes and Firmicutes further signals a substantial difference in the oxygen conditions between fractures and matrix sediments, and supports a conversion from primarily aerobic to anaerobic conditions between 230 and 325 cmbgs. Within the Bacteroidetes, which increase with depth in accordance with a previous study (Pereira et al. 2017), Cytophagaceae and Chitinophagaceae are families found primarily in the upper fractures. The few isolated strains from these families have been characterised as aerobic (Chaudhary, Dahal and Kim 2018; Li et al. 2018). In contrast, the anaerobic Bacteriodaceae (Satoh et al. 2002) become more abundant with depth. For the Firmicutes, facultative anaerobic Bacilli were dominant in the upper samples, but replaced by obligate anaerobic Clostridia in deeper matrix samples. For comparison, Pereira et al. (2017) found Bacilli to be the dominant class of Firmicutes at all depths down to 700–800 cmbgs in a sandy soil with a higher porosity than the current clayey till. Clostridia contain sulphate reducers that are more pronounced in saturated soils (Hansel et al. 2008), which may explain why these were detected in the current deep sediments while they were not detected by Pereira et al. (2017). Additionally, some Deltaproteobacteria are known sulphur reducers (Hansel et al. 2008), corresponding well with the high relative abundance of this class at 465 cmbgs in both fractures and matrix sediments and the detection of high relative abundances of dsrAB genes at the same depth. The lack of dsrAB genes in the upper part of the soil has also been seen in a previous study of unsaturated soil classified as Inceptisol (Hansel et al. 2008).

The relative abundances of genes involved in aerobic as well as anaerobic N cycling were also used as proxies for oxygen availability. Furthermore, the abundances of amoA, nirK and nirS genes have previously been linked to the corresponding process rates of nitrification and denitrification (Rocca et al. 2015). An assumed influence of oxygen was observed for amoA genes involved in ammonia oxidation, showing the highest relative abundance in the upper layers and occurring in deeper preferential flow path samples than matrix sediment samples. Although oxygen levels were not measured, visual inspection showed precipitated iron and manganese oxides on the deeper fracture walls (Fig. 1), confirming their aerobic status (Beven and Germann 1982).

In contrast, the nirK and nosZ genes involved in denitrification became more prevalent with depth. Denitrifiers containing nirK are highly diverse, and it is difficult to link gene abundance with 16S rRNA gene profiles (Zumft 1997). However, the detection of nosZ genes could potentially be linked to alphaproteobacterial Rhodospirillales (Fig. S4) (Harter et al. 2017), or to taxa within the Betaproteobacteria and Gammaproteobacteria as these groups are targeted by the primers used for the nosZ1 gene (Henry et al. 2006). The abundance of nitrite reductases in deep clay tills has not previously been studied, but a higher relative abundance of nirK and nirS genes in relation to 16S rRNA genes was observed in sediments below rice paddies at 100 cmbgs (Wang et al. 2017), probably due to the prevailing anaerobic conditions and a higher input of nitrate. In general the nirS genes had a higher relative abundance than the nirK genes, which is in agreement with the findings in river sediments (Huang et al. 2011), topsoil (Tsiknia et al. 2015) and subsoil (Zhu et al. 2018). The current distribution of amoA, nirK and nosZ genes suggests the co-occurrence of nitrifiers and denitrifiers, which is indicative of environments with fluctuating aerobic and anaerobic conditions. The gene distributions point to a major shift from aerobic to anaerobic N cycling at ~325 cmbgs in the matrix in agreement with the sediment changing colour from reddish brown (oxidised) to bluish grey (reduced) (Fig. 1). However, this transition occurred at a greater depth in the preferential flow paths.

Interestingly, the distribution of the order Gallionallales within the Betaproteobacteria showed a high abundance in the deep matrix sediment at 465 cmbgs, while this group was almost absent from the adjacent fracture. Gallionallales oxidise iron under microaerophilic conditions (Emerson, Fleming and Mcbeth 2010) and their occurrence indicates a specific redox boundary layer required for bacterial iron oxidation at this depth.

In the matrix sediment, heterotrophic Planctomycetes primarily appeared down to 95 cmbgs. By comparison, it was abundant in the preferential flow paths at all depths, even in the sand lens at 575 cmbgs. This phylum is recognised for its ability to aerobically degrade plant debris including saccharides and heteropolysaccharides (Erbilgin, McDonald and Kerfeld 2014; Wang et al. 2015), and these bacteria probably benefit from the input of organic matter provided by the preferential flow paths. Even other taxa, such as Chloroflexi, with abilities to degrade complex plant materials, were enriched in preferential flow paths (Seuradge, Oelbermann and Neufeld 2017). This suggests that the supply of organic matter in the preferential flow paths have an impact on the bacterial communities, and also the total abundance of bacteria.

According to the NMDS analysis of the bacterial communities, communities from the matrix sediments and corresponding preferential flow paths generally did not cluster, reflecting the very complex nature of the deep sediment environment and their bacterial communities. However, there were two interesting observations of clustering. First, communities from the plough layer and the underlying biopores clustered together and were both characterised by a high Shannon diversity index. The close relationship between these two communities, also supported by the nestedness analysis, can be attributed to the strong effect earthworms have on microbial communities (Banfield et al. 2017). A high diversity of communities from linings of biopores (drilosphere) has also been observed for a German clayey till (Uksa et al. 2015). However, the biopores in the current study were sampled from a greater depth (70–120 cmbgs) than previous studies (Stromberger, Keith and Schmidt 2012; Uksa et al. 2015), expanding the influence of earthworm activity to bacteria in deeper sediments. Earthworm burrows were observed down to 1.5 mbgs, and the role of biopores for bacterial communities probably becomes more important with depth because of the limited supply of available C (Hoang et al. 2016).

Second, we observed a clustering of communities in the sand lens at 575 cmbgs and in the matrix sediment at 95 cmbgs. In contrast to communities in biopores, communities in tectonic fractures and sand lenses have only been studied to a very limited extent. Surprisingly, these communities had a high abundance of bacteria and diversity, pointing to a bacterial hotspot environment at a relatively great depth. The clustering of the sand lens samples with top soil communities for both 16S rRNA genes (Fig. 3) and functional genes (Fig. 5), and the detection of archaeal amoA and nirK genes accentuate the sand lens as a soil compartment with a high potential for both aerobic and anaerobic nutrient cycling, which deserves further investigation in the future.

It was shown that communities in some preferential flow paths had more indicator OTUs from the plough layer than corresponding matrix sediments. Furthermore, the applied nestedness analysis determined that bacteria from the plough layer and biopores influenced the bacterial communities in preferential flow paths more than matrix communities. As discussed above, the resemblance of communities in the plough layer and biopores could be explained by the impact of the earthworms (Stromberger, Keith and Schmidt 2012), but we also speculated that water-driven transport of bacteria could play a role in the composition of communities in deeper preferential flow paths. As previously reported, root-associated bacteria in topsoil (Dibbern et al. 2014) and E. coli. (Arnaud et al. 2015) are transported to the deeper layers of the soil and groundwater aquifers. Although active bacterial transport was not measured, a greater relative abundance of Rhizobiales was found in the biopores, fractures at 230 cmbgs and the sand lens than in the corresponding matrix sediments. Roots may prefer to grow within fractures (Jørgensen et al. 2002), and it is possible that root-associated Rhizobiales are transported to greater depths via roots. Our findings strongly suggest that subsoil preferential flow path communities are influenced by down-passing bacteria from the plough layer communities above. The detection of many indicator OTUs from the preferential flow paths and plough layer in the sand lens, combined with community composition and increased diversity and abundance, identify sand lenses as potential aerobic hotspots in the clayey till.

CONCLUSIONS

The current study reinforces the importance of including preferential flow paths of clay till in future studies since they represent an almost unexplored niche of the subsoil ecosystem. It also highlights communities in biopores, fractures and sand lenses with higher diversity and a different composition than those in matrix sediments. The increased abundances of aerobic taxa indicate that aerobic conditions extend to a greater depth in the fractures than in the adjacent matrix sediment, while greater relative abundances of plant-material degraders point to the influence of transported organic matter on community composition in the flow paths. In addition to the effects of abiotic conditions on bacterial communities in flow paths, indicator OTU and nestedness analyses suggest an impact of vertical transport of bacteria.

Furthermore, the study points to the significant effects of depth on bacterial community compositions of both preferential flow paths and matrix sediments. While aerobic bacteria and genes involved in aerobic nutrient cycling dominated the upper sediments in both matrix and preferential flow paths, anaerobic bacteria and genes involved in anaerobic nutrient cycling became dominant in deeper sediments. Nevertheless, preferential flow path habitats as deep as 575 cmbgs harboured potential for both aerobic and anaerobic nutrient cycling. Hence, it is suggested that preferential flow paths should be included when evaluating microbial diversity and modelling nutrient turnover in clayey tills.

Click here for additional data file.