Response of Archaeal and Bacterial Soil Communities to Changes Associated with Outdoor Cattle Overwintering

Archaea and bacteria are important drivers for nutrient transformations in soils and catalyse the production and consumption of important greenhouse gases. In this study, we investigate changes in archaeal and bacterial communities of four Czech grassland soils affected by outdoor cattle husbandry. Two show short-term (3 years; STI) and long-term impact (17 years; LTI), one is regenerating from cattle impact (REG) and a control is unaffected by cattle (CON). Cattle manure (CMN), the source of allochthonous microbes, was collected from the same area. We used pyrosequencing of 16S rRNA genes to assess the composition of archaeal and bacterial communities in each soil type and CMN. Both short- and long- term cattle impact negatively altered archaeal and bacterial diversity, leading to increase of homogenization of microbial communities in overwintering soils over time. Moreover, strong shifts in the prokaryotic communities were observed in response to cattle overwintering, with the greatest impact on archaea. Oligotrophic and acidophilic microorganisms (e.g. Thaumarchaeota, Acidobacteria, and α-Proteobacteria) dominated in CON and expressed strong negative response to increased pH, total C and N. Whereas copiotrophic and alkalophilic microbes (e.g. methanogenic Euryarchaeota, Firmicutes, Chloroflexi, Actinobacteria, and Bacteroidetes) were common in LTI showing opposite trends. Crenarchaeota were also found in LTI, though their trophic interactions remain cryptic. Firmicutes, Bacteroidetes, Methanobacteriaceae, and Methanomicrobiaceae indicated the introduction and establishment of faecal microbes into the impacted soils, while Chloroflexi and Methanosarcinaceae suggested increased abundance of soil-borne microbes under altered environmental conditions. The observed changes in prokaryotic community composition may have driven corresponding changes in soil functioning.


Introduction
Upland grassland soils are usually well aerated, characterized by high levels of dissolved organic carbon (C) and acidic pH. Microbial growth in these systems is mainly limited by the availability of nitrogen (N) and phosphorous (P). These conditions favor the Acidobacteria, Proteobacteria, Actinobacteria, and Firmicutes [1][2][3], and archaea such as Thaumarchaeota [4]. Grassland soil conditions undergo significant changes when used for intensive grazing and husbandry. Such sites receive large inputs of N in the form of urine and manure [5], which also increases soil pH [6,7]. In addition, soils become compacted, depending on the number of animals per hectare, affecting oxygen availability and soil redox conditions. Intensive grazing may cause shifts in microbial communities associated with changes in soil physicochemical properties [8][9][10][11][12] or due to the introduction of faecal microbiota [9]. Ultimately, processes, such as denitrification and methanogenesis, increase dramatically, particularly when overall microbial activity is high at sites with high animal traffic [5,8,[13][14][15].
Outdoor winter husbandry has become more popular in Europe over recent decades, partly as a result of the introduction of open cattle sheds. Compared to summer grazing, animal traffic is much higher over a smaller area, often situated close to the cattle shed. In addition, winter soils undergo frequent freezing and thawing cycles and overall soil water contents is usually higher, while vegetation cover is scarce or non-existent. Under these conditions, it is expected that soil microbe responses would be amplified and would lead to large community shifts.
The strong impact of outdoor winter husbandry on microbial groups that drive redox processes, such as the denitrifiers [13] and methanogens [8], affect the emissions rates of greenhouse gases such as nitrous oxide (N 2 O) and methane (CH 4 ), with levels increasing strongly over summer grazing pastures [5,8,13]. This raised the question whether the overall microbial diversity responds to the form of soil management. Indeed, in a recent study using phospholipid fatty acids (PLFA) as microbial biomarkers [10], we were able to show that overall soil microbial biomass increased as a consequence of winter outdoor husbandry, with all groups investigated (i.e. bacteria, fungi and archaea) affected, though to differing extents. This was confirmed through fingerprinting based methods on selected ribosomal genes for archaeal, bacterial, and fungal communities [9]. However, a comprehensive evaluation of the changes on microbial community composition was not possible using this approach.
In the present study we aimed to investigate the composition of archaeal and bacterial communities using pyrosequencing of 16S rRNA amplicons, thereby identifying major responders to the environmental changes caused by cattle overwintering. Therefore, differences in soil properties were tested and correlated with changes in soil microbial community composition. We also aimed to evaluate the resistance and resilience of soil prokaryotic communities by looking at areas submitted to short-and long-term cattle impact and a naturally regenerating area. In addition, we assessed the archaeal and bacterial composition of the cattle manure in order to determine if microbes entering the soil via manure may putatively survive in this environment. Based on data obtained in previous study [9] we expect to observe different response pattern for archaeal and bacterial communities. We hypothesize that cattle impact will alter the soil microbial community composition with negative effects on diversity relative to the control. Moreover, we expect to observe significant changes in relative abundance of taxa sensitive to cattle overwintering associated disturbance (e.g. oligotrophic vs copiotrophic taxa, aerobic vs anaerobic microbes, introduction of faecal microbes) and specific functional guilds (e.g. methanogens vs ammonia-oxidizers).

Site description
Soil samples were collected from a cattle farm near the village Borová in the Czech Republic (48°54' 51" N, 14°14' 51" E). No specific permissions were required for these locations and activities. The farm is private land; the owner is Mr. Kamír, who should be contacted before future access. During the field work, there was no contact with animals as we collected only cattle manure deposited on the soil surface. Therefore no approval was obtained. The site is characterized by a mean annual temperature of 7°C and precipitation of 650 mm. The predominant soil type is sandy loam, classified as Haplic Phaeozem (arenic; WRB system). From October to May each year since 1995, approximately 90 cows have been moved to an open cattle-shed connected with a 4 ha overwintering area (during the summer vegetation period, the animals graze on pastures located elsewhere). In general, the impact of the cattle is greatest close to the cattle shed. In May 2011, we collected soil samples from four locations on the farm representing different levels of impact. The first location, representing soil under long-term impact (LTI), was located within a 25 m radius of the cattle shed and has been used continuously for overwintering since 1995. The second, representing soil regenerating from moderate impact (REG), was located in a pasture close to the cattle shed that had been used for overwintering each year between 1995 and 2008. The third location, representing soil under short-term impact (STI), was situated in a pasture more than 50 m away from the cattle shed that has been used for overwintering since 2008. A control sample (CON) was located outside of the pastures and was dominated by perennial grasses and clovers (S1 Fig).
The upper horizon of LTI soil is actually a mixture of soil and cattle manure. It is a dark brown, rich in organic matter and vegetation is heavily damaged. Porosity is relatively high and the topsoil (0-20 cm) bulk density is lower than that in the CON soil due to the high input of organic matter. In contrast, the lower soil horizon is more compact due to long-term disturbance and trampling. STI soil reflected similar changes regarding the vegetation cover and mixing of soil upper layer with cattle manure, but less pronounced than in LTI. On the other hand, REG soil does not show any visible differences from CON since vegetation is recovered.

Soil and cattle manure sampling and analyses
Five independent composite samples were taken in May 2011, each consisting of seven subsamples randomly collected over an area of 1 m 2 . Single subsamples were taken from the upper 20 cm using a 5 cm diameter auger. In order to reduce heterogeneity within the composite samples, subsamples were well mixed and sieved (5 mm mesh) on site. Replicates at each location were taken at least five meters away from the previous sample. Therefore, plots were non-randomized across the farm area and our replicates could be considered pseudo-replicates [16]. On the other hand, given the uniformity of the landscape, we may assume that the four plots had similar soil properties prior to the establishment of the overwintering pasture 15 years ago. In other words, it is reasonable to assume that the effects of treatments on microbial parameters were more determined more by cattle impact than initial differences in soil properties across plots. Five composite CMN samples were taken at the same time and analyzed separately, each comprised from subsamples taken from ten individual cattle manure deposits on the field. To reduce heterogeneity subsamples were homogenized. Each soil and CMN composite sample was divided in two parts; one was transferred into cryovials and kept at dry ice and finally stored at -80°C until DNA extraction, whereas the second one was stored at 4°C for analysis of physicochemical properties.
Major physical and chemical properties of the four soil types and their correlations are outlined in Table 1 (raw data are available in S1 Table) and S2 Table. Soil organic carbon (Corg) was determined by wet oxidation with acid dichromate; total nitrogen (Ntot) was measured by Kjeldahl digestion [17]. Gravimetric moisture content (GWC) was determined after drying at 105°C for 48 h. pH was determined using a pH meter in 1:2.5 (w/w) soil/CMN: destilled water suspension. Cation exchange capacity (CEC), P, K, Ca, and Mg concentrations were measured as described by Zbíral et al. [18]. The same techniques were applied to measure physical and chemical properties of CMN samples.

DNA extraction
DNA was extracted from 0.5 g of fresh soil and CMN samples using the FAST Spin kit for Soil (MP Biolabs, USA). The quality and quantity of DNA was checked through gel electrophoresis and use of the PicoGreen assay (Invitrogen, USA).
All samples were amplified in triplicate, pooled, and purified using the QIAquick PCR Purification kit (Qiagen). Amplicon quality was checked using agarose (2%) gel electrophoresis and quantified using the Quant-iT dsDNA BR assay kit (Invitrogen, USA). The amplicons were pooled equimolar to create two forward and two reverse oriented libraries for both Bacteria and Archaea.

Pyrosequencing
Four single-stranded libraries with different sample-specific adaptors were cleaned up from smaller fragments within each sample using Agencourt AMPure beads (Agencourt Bioscience Corporation, MA, USA); quality was evaluated using a Bioanalyzer 2100 (Agilent Technologies, Germany) and DNA 1000 LabChip software. Single-stranded DNA libraries were generated using the GS FLX Standard DNA Library Preparation Kit (Roche). Uniquely tagged, pooled DNA samples were immobilized onto DNA capture beads, amplified through emulsion-based clonal amplification (emPCR), and sequenced using Titanium reagents and procedures in four regions of a PicoTiterPlate on a 454 Genome Sequencer FLX system (average

Pyrosequencing data analysis
Initial signal processing and quality filtering of pyrosequencing reads was performed using the automatic amplicon pipeline of the GS Run Processor (Roche) in order to remove failed and low-quality reads from the raw data. Sequencing errors and chimeras were removed using the programs Amplicon noise and Perseus according to Quince et al. [23]. Forward and reverse DNA sequences with an exact match over at least 100 bp were assembled via the custom C program using exact pairwise Needleman-Wunsch alignments [24]. Sequences were clustered into Operational Taxonomic Units (OTUs; defined as a group of sequences sharing 97% nucleotide sequence identity) using UCLUST software [25]. Singletons were removed from the dataset. Taxonomic classification and assignment of individual OTUs was performed using CREST LCAClassifier against the SilvaMod SSU rRNA reference database ( [26], http://apps.cbu.uib. no/crest). The similarity cutoff for assignment has been set as follows: for the genus, family, order, class and phylum ranks the respective cut-offs are 97%, 95%, 90%, 85% and 80% identity, respectively [26].

Nucleotide sequence accession number
The 16S rRNA gene sequences derived from pyrosequencing have been deposited in the NCBI Sequence Read Archive under accession number SRP041238 (BioProject No. PRJNA240697).

Statistical analysis
All diversity analysis were calculated (based on 97% nucleotide sequence identity) using the vegan statistical software package for R (v.2.15.1; [27]). We determined alpha diversity indices, estimated using (i) species richness, calculated as Margalef index, (ii) the Shannon-Weaver index (H') and (iii) the Pielou evenness index (j). Sørensen index (S), a percent similarity measure, for which larger numbers indicate greater similarity, was used for the computation of the beta diversity. Effect significance was determined through nonparametric Kruskal-Wallis analysis of variance (ANOVA) followed by the post hoc Tukey HSD test. Significance of correlations between physicochemical properties and diversity indices, and between soil edaphic factors and taxon abundance, was tested using Pearson's correlation coefficient using the vegan statistical software package for R (P >0.05). Heatmaps for all OTUs were calculated using the gplots package for R [28]. Multivariate statistical analysis was used in order to explain variation in the data and to test the significance of cattle overwintering on both soil chemistry and archaeal and bacterial community composition. Abundance data were transformed using log (x+1), centered, and standardized by total (absolute values converted to relative values) prior to construction of a Bray-Curtis similarity matrix. Physico-chemical properties (except pH) were also log (x+1) transformed. Detrended canonical correspondence analysis (DCCA) was used to determine the length of gradient along the first ordination axis in order to select the appropriate ordination method for the data. Distance-based Redundancy Analysis (db-RDA) was performed to assess the relationship between known environmental variables and variation in the multivariate data. The significance of the relationship was tested using the Monte Carlo permutation test [29], with 499 or 999 unrestricted permutations and manual forward selections used for the db-RDA. Because of high collinearity between explanatory variables, the variation inflation factor (VIF), which expresses the extent of multiple correlations with other predictors, was measured. Variables with VIF<20 were added to the final db-RDA analysis. All analyses were performed using the CANOCO 5 software package [29]. To study effects of grazing intensities on community composition and if communities differed between studied soils we conducted a PerMANOVA based on Bray-Curtis dissimilarity matrices and pair-wise comparisons [30] using the PRIMER-6 package [31]. Analysis of similarity (ANOSIM) as well as similarity percentage (SIMPER) analysis was used to determine the level of similarity between the samples and to identify the taxa that were mainly responsible for the differences observed between soil samples with different grazing intensities.

Archaeal and bacterial community richness, diversity, and evenness
A total of 1,092,472 reads were obtained with an average length of 496 bases. Half of these sequences (527,820) passed the quality filters and were subsequently used for further analysis, 247,832 derived from bacterial amplicons and 279,988 from archaeal amplicons. The number of classified sequences per sample ranged from 6,458 to 9,891 for Bacteria and from 6,096 to 11,304 for Archaea, respectively. Following assembly of forward and reverse reads in each sample, we obtained 4,415 to 5,843 sequences for Bacteria and from 7,990 to 10,471 sequences for Archaea. As the number of sequences in different samples did not differ significantly ( Table 2), no further subsampling was undertaken.
OTU saturations at the 97% nucleotide sequence identity level were obtained for bacterial and archaeal data sets when 4,000 and 1,000 sequences were sampled, respectively (S2 Fig). Table 2 shows the diversity analyses carried at OTU 97% for both data sets. For archaeal communities the lowest α-diversity values (H', R, and j) were obtained for STI. Archaeal α-diversity Table 2. Richness estimates, diversity and evenness indices of soil (LTI = long-term impact; REG = regenerating soil; STI = short-term impact; CON = control) and cattle manure (CMN) archaeal and bacterial communities based on OTU clustering at 97% nucleotide sequence identity. was negatively correlated with organic C, while richness was additionally negatively related to pH, CEC, and total N ( Table 3). Bacterial community diversity and richness did not significantly differ among the investigated soils. However, bacterial evenness was significantly lowest in LTI soil ( Table 2). Strong correlations were observed between bacterial evenness and CEC, organic C, pH, and total N as well as between bacterial α-diversity and CEC, and organic C ( Table 3). Furthermore, β-diversity (Sørensen index; Table 2) indicated greater similarities among bacterial communities than among archaeal communities. This was confirmed by the heatmaps (S3 Fig), which showed differences between samples based on relative abundance of individual OTUs in archaeal and bacterial communities. However, in both cases, the lowest βdiversity indices were obtained in the LTI soil.

Archaeal community composition
We obtained 55 archaeal OTUs, all of which were affiliated to one of three phyla, i.e. Thaumarchaeota, Euryarchaeota, and Crenarchaeota. None of the sequences remained unclassified at phylum level (cutoff 80% sequence similarity [26]  high proportion of Methanobrevibacter as well as Methanocorpusculaceae and vadinCA11 gut group ( Table 4).
The remarkable difference between CMN and soil archaeal communities was verified by db-RDA ( Figure A in S4 Fig). In accordance to ANOSIM, db-RDA ordination of soil communities based on OTU-level showed clustering of soils according to pre-defined groups (Fig 1A). The explanatory variables accounted for 89.5% of variation (total variation 3.26, Pseudo-F = 32.1, P = 0.002). Iterative forward selection resulted in Ca concentration and grazing   (Fig 1A).
These findings were later confirmed by co-occurrence analysis ( Figure A in S5 Fig). Methanocorpusculum, Methanosphaera, Methanobrevibacter, and TMEG-2 OTUs overlapped in CMN and LTI soils, while only Methanosphaera, and Methanobrevibacter were shared by STI soils and CMN. Candidatus Nitrososphaera and members of SCG I.1b and Group I.1c were only detected in REG and CON soils. By comparing CON and STI soil archaeal communities we aimed to identify major responders to short term overwintering grazing. SIMPER analysis (dissimilarity contribution > 3%) showed that members of SCG group I.1b (OTUs 10, 13, 14) and Group I.1c (OTUs 15, 16) significantly declined and Methanobrevibacter (OTUs 2, 3) significantly increased in STI soil. On the other hand, looking at dissimilarities between REG and LTI identified archaeal responders associated with regeneration of grassland soil from long term cattle impact. Two SCG I.1b members (OTUs 0 and 1) increased and two MSC members and Methanosarcina significantly decreased in REG soil in comparison to LTI (OTUs 44, 51, and 7, respectively).

Bacterial community composition
Overall, 2,591 bacterial OTUs from 25 phyla were identified. Of these, 24 phyla were detected in soil samples and 18 in CMN ( Table 4). Between 9 and 14% of sequences remained unclassified as they were below the 80% classification threshold ( [26]; see Materials and Methods). CON and REG soils were characterized by a dominance of Acidobacteria and Proteobacteria, together representing 81.8% of all CON and 72.8% of all REG sequences ( Table 4). Next, Actinobacteria and Bacilli (Firmicutes) together represented 9.8% of the CON community and 10.6% of the REG community. Acidobacteria negatively responded to cattle impact and followed the trend CON REG >STI > LTI, with the relative abundance significantly lowered in STI and LTI compared to CON (P < 0.05, Table 4). Similar trends we observed for Planctomycetes, Proteobacteria (α-and βclasses), and Verrucomicrobia. The Firmicutes followed the opposite trend, being 2 to 3 times higher in STI and LTI soils, respectively, in comparison to CON. The largest positive response to long term cattle overwintering was observed for Clostridia, which reached the abundance of 12.4% in LTI, while in CON soil were rare (0.2%). Moreover, Bacteroidetes, Deinococcus-Thermus, Fibrobacteres, Bacilli, Gemmatimonadetes, Chloroflexi and γ-Proteobacteria were also more abundant in cattle impacted soils. STI was enriched by Firmicutes (both Clostridia and Bacilli), Gemmatimonadetes, Chloroflexi, and β-Proteobacteria in comparison to CON, whereas Acidobacteria, Planctomycetes, αand γ-Proteobacteria, and Verrucomicrobia were depleted. REG bacterial community showed recovery of Acidobacteria, Planctomycetes, and α-Proteobacteria in comparison to LTI, while Bacteroidetes, Firmicutes, Chloroflexi, Deinococcus-Thermus, Fibrobacteres, Verrucomicrobia, and γ-Proteobacteria decreased. Bacterial communities in four studied soils significantly differed from each other albeit the strength of test was lower than in case of archaeal counterparts (PERMANOVA, Pseudo-F = 8.53, p = 0.001, 996 permutations; ANOSIM, Global R = 1, P = 0.001). SIMPER pairwise test showed the clustering of the soil bacterial communities according to grazing intensities (respective average dissimilarities (%): CON-REG = 58.91, CON-STI = 71.13, CON-LTI = 73.69, REG-STI = 41.62, REG-LTI = 54.16, STI-LTI = 59.36). SIMPER analysis showed large dissimilarity of between CMN and soil bacterial communities, expressing average dissimilarities in the range of 97.4-99.75%. This was confirmed by db-RDA ( Figure B in S4 Fig).
In db-RDA, explanatory variables accounted for 57.9% of variation (total variation 4.64, Pseudo-F = 3.3, p = 0.002, Fig 1B). The PC1 axis was characterized by increases in soil water content, pH and CEC, organic C, total N, P, K, magnesium (Mg), and Ca concentrations, all of which favoured the LTI soil bacterial community. Iterative forward selection resulted in grazing intensities and pH as significant explanatory variables, contributing by 38.7% and 16.6% to the model and explained 23.0% and 6.2% of variability, respectively. The PC2 axis discriminated REG and STI bacterial communities from that of CON and LTI. Differences between STI and REG were marginal and the lowest among the soils studied, supporting the SIMPER results (average dissimilarity < 50%). These results were additionally confirmed by co-occurrence analysis ( Figure B in S5 Fig). The 30 OTUs displaying highest ordination correlation were plotted into the diagram (Fig 1B). Rhizomicrobium and three other OTUs of the Acidobacteriaceae group 6 (Ac 6) and β-Proteobacteria were negatively correlated with the PC1 axis, while most others were clustered at the opposite pole, indicating dominance in LTI. Diverse anaerobic bacteria were identified among them, including Lachnospiraceae, Propionibacteriaceae, Chloroflexi (including Anaerolineae), Thauera, Proteiniclastisum, and the facultative anaerobe Caldilinea. Some of the LTI-correlated OTUs, such as Caldilinea, Devosia, Truepera, and Dietzia, were typical of alkaline habitats. Differences in relative frequencies of the most abundant bacterial genera are shown in S3 Table. SIMPER analysis showed that 812 OTUs (data not shown) contributed to dissimilarity between CON and STI bacterial communities. Three major OTUs (dissimilarity contribution > 0.5%) were assigned as Planococcaceae, Comamonadaceae, and Pseudoxanthomonas (OTUs 3, 2, and 0, respectively). Similarly, 671 OTUs contributed to difference in the bacterial community structures between REG and LTI soils. Among them, Ac 6 members (OTUs 5 and 14) and Pedomicrobium (OTU310) significantly increased and different Firmicutes (Lactobacillales-OTU7, Proteiniclasticum-OTUs 17 and 350, Anaerolineaceae-OTU2242) and Actinobacteria (Corynebacterium, OTU20) members decreased in REG soil in comparison to LTI.

Archaeal and bacterial community changes at higher taxonomic levels
At higher taxonomic level (archaeal phyla and classes, bacterial phyla) the community shifts visualized by db-RDA ordinations were characterized by the similar trends as for respective OTU-levels. In Archaea, the explanatory variables (grazing intensity and Ca concentration) together accounted for 94.5% of variation (total variation 1.29, Pseudo-F = 64, p = 0.002, Figure A in S6 Fig). The PC1 axis separated LTI and REG soils from STI and CON soils, with Ca concentration showing high correlations; thus indicating that the LTI soil community had increased concentrations of Ca and other nutrients, these being highly correlated with each other (P < 0.05, S2 Table). Relative abundance of most of archaeal lineages was strongly correlated also to soil pH (S4 Table). The largest positive response were found for Methanobacteriaceae, while the strongest negative correlations were found for the Marine Group II, Group I.1c, and SAGMCG-1 lineages. In contrast to higher taxonomic levels, most of methanogenic genera responded strongly to organic C and CEC (S4 Table).
Db-RDA based on bacterial phyla level, showed higher difference between STI and REG than respective analysis based on OTU abundance. In contrast, soil pH did not contribute to the model significantly, thus grazing intensity was the only explanatory variable. In addition, it indicated that Acidobacteria, Elusimicrobia, Planctomycetes, Verrucomicrobia, and Candidate divisions OP11 and WS3 were strongly correlated with CON soil (total variation = 0.43, explanatory variables accounted for 53.8%, Pseudo-F = 6.2, p = 0.002; Figure B in S6 Fig). Bacterial indicator groups for LTI soil included Bacteriodetes, Cyanobacteria, the Deinococcus-Thermus Group, Fibrobacteres, Lentisphaerae, Synergistetes, Tenericutes, and Candidate Division OD1. REG soil was correlated with Nitrospirae, while STI soil was associated with Gemmatimonadetes, and Erysipelotrichi ( Figure B in S6 Fig). Bacterial group relative abundance was significantly correlated with soil pH and other physicochemical soil properties (i.e. total N, organic C, and CEC) at a range of taxonomic levels (phylum, class, or genus) (Spearman's rank correlation, P > 0.05; S4 Table). Abundance of Acidobacteria, Proteobacteria (mainly α-Proteobacteria, including Bradyrhizobium, and Bryobacter and Sporosarcina), and Planctomycetes decreased with increasing pH and physicochemical factors. In contrast, Firmicutes (both Clostridia and Bacilli), Chloroflexi, Actinobacteria, Elusimicrobia, Deinococcus-Thermus, Candidate phylum OD1, and gg. Anaerolinea, Caldilinea, Devosia, Dokdonella, Peptostreptococcaceae Insertae Sedis, Proteiniclasticum, Tetrasphaera, and Trichococcus are all increased with pH.

Different responses of soil archaea and bacteria to cattle overwintering
In this study, we reported the significant negative response of soil archaeal community, altering both diversity and composition, to cattle outdoor husbandry. The effect on soil bacterial αdiversity was of less extent, while shift in bacterial community composition was significant. These results corroborate our previous findings observed by DGGE (denaturing gradient gel electrophoresis) fingerprinting of 16S rRNA gene amplicons [9], indicating higher resilience of soil Bacteria to cattle overwintering induced changes in comparison to soil Archaea. Significant loss of beta diversity found in long-term impacted soil suggests that cattle overwintering leads to the homogenization of both microbial communities over time. While archaeal diversity was generally low in our study, we observed significantly higher diversity than was reported in a previous study [32], in which clone library and sequencing were adopted. Auguet et al. [4] were the first to report unexpectedly low soil archaeal diversity and richness while searching for global ecological patterns in Archaea. Despite calculating diversity at the lineage level (a higher taxonomic level than OTU), their results showed few archaeal lineages present in soils, and few of these as dominant.

Shifts in archaeal communities
Shifts in archaeal community composition were influenced by the intensity of cattle overwintering and by increased concentration of Ca (identified by manual forward selection, correlated to other factors, see S2 Table). Thaumarchaeota (formerly described as Crenarchaeota) dominated archaeal community of CON soil. This confirms previous study [4] that identified Thaumarchaeota I.1b and I.1c as dominant indicator lineages for the soil environment, other lineages showing only moderate (Methanosarcinales, Methanomicrobiales, Thermoplasmatales, and Halobacteriales) or very low (Thaumarchaeota I.1a, C2, Methanobacteriales) abundance. Both Thaumarchaeota SCG-I.1b clade, dominating in CON, and SAGMCG-1clade are assumed to contribute to ammonia oxidation in well aerated soils, especially at neutral and acidic pH. Relatives of Groups I.1b and I.1c have also been identified in natural and managed grasslands [32] and in other soils with a broad pH range [24]. In upland pastures, increasing pH, associated with deposition of animal urine and manure, is strongly linked with shifts in the archaeal community [9,12,33,34]. In our study, Thaumarcheota differed in response to overwintering grazing. In particular, members of SCG-I.1b remained dominant in the community only in REG and STI soils, while size of Group I.1c and SAGMCG-1 declined in all cattle impacted soils. This is in line with the results of Gubry-Rangin et al. [24], who recognized Group I.1b as the dominant Ammonia Oxidizing Archaea in all soils along a broad pH gradient. Group I.1c and SAGMCG-1 are associated with a range of acidic environments [35][36][37][38][39] and showed negative correlation to soil pH, indicating a possible inability to adapt to increased pH. Both groups showed the highest sensitivity to environmental changes associated with cattle outdoor husbandry, being absent in other soils. On the other hand, SCG I.1b lineage showed resilience at least to moderate cattle impact represented here by STI, while become sensitive under long-term practice.
The proportion of euryarchaeal sequences, mainly methanogenic, increased significantly with intensity and duration of cattle impact, reflecting the significant changes in soil physicochemical properties. Euryarchaeal sequences in CON soil were associated with the non-methanogenic Thermoplasmata Marine group II. In previous studies, methanogens have been detected in CON, but at abundance levels around the limit of detection (ca. 10 5 mcrA gene copies per g dw soil, [8,15,40]). Not surprisingly, consumption rather than emission of methane was detected in CON soil [8]. The archaeal community of LTI soil differed significantly from that of other soils, being characterized by reduced Thaumarchaeota, an increase in MCG and a dominance of diverse methanogenic members. MCG have been detected in a range of marine and continental habitats, and are widely distributed in subsurface anoxic and semi-anoxic habitats [41,42]. The ecological role of MCG remains unclear; however, in the light of recent genomic insights, it is thought these archaea may contribute significantly to degradation of recalcitrant organic matter, including aromatic compounds [43]. The fact, that that MCG are indicative for long term cattle overwintering, may implicate that they favour nutrient rich, alkaline soils with increased anaerobic niches. The abundance of methanogens, Methanosarcinaceae (41%) and Methanomicrobiaceae (1.5%) in particular, increased in LTI soil as a result of altered environmental conditions. This corroborates our previous findings of significant methanogenic activity, dominated by Methanosarcina and uncultured rumen archaea, in LTI soil [8,11]. Methanosarcina is versatile, following all three methanogenic pathways, and favor high pH (optimum 8.0 [44]) and acetate concentrations. Thus, Methanosarcina should be favored in soils with a high input of organic material in form of cattle slurry, where higher production of acetate [45], acetic acid and methylamines [46][47][48] may occur. In REG soil, the methanogenic community was similar to that of LTI soil, but significantly less abundant (<5%). On the other hand, STI soil was enriched by Methanobacteriaceae, likely introduced with cattle manure, being the most abundant methanogen group in CMN (up to 76%; Table  4). Since Methanobacteriaceae sensitively reacts to changes in cattle overwintering practice, their increased abundance could act as an indicator of altered soil properties associated with this management. CMN samples comprised Methanobacteriaceae, Methanocorpusculaceae, and the vadinCA11 Gut Group (vertebrate intestinal cluster), which corroborates previous findings describing the archaeal community of cattle rumen or manure [48,49].While the transfer of methanogens from cattle intestinal tracts into soil has been reported previously [8,50], members of the Methanocorpusculaceae and vadinCA11 (Methanoplasmatales) lineages were not detected in our soils, despite being highly abundant in CMN. This finding corroborated previous studies suggesting that both groups do not survive outside the intestines. Differences between the methanogenic community of CMN and cattle-impacted soils (mainly LTI) reflect changes in nutrient supply and the presence of available electron donors for methanogenesis. While rumen methanogens (e.g. Methanocorpusculaceae and Methanobacteriaceae, [49]) mostly utilize hydrogen, in soils other sources can also be utilized effectively (e.g. the formate, acetate, and methyl groups). Moreover, typical rumen methanogens are often strongly associated with anaerobic protozoa conducting interspecies hydrogen transfer (endosymbionts; [51,52]) and/or they are attached to the rumen epithelium [53] and, therefore, may be tightly adapted to a specific niche in the cattle intestine.

Shifts in bacterial communities
This study confirmed that the intensity of cattle husbandry have a significant effect on soil bacterial community composition, as documented by changes in the relative abundance of bacterial phyla, classes, and OTUs in the different soil types. Grazing intensities and soil pH were recognized as factors explaining the most of the variability in soil bacterial community composition. Most studies on the composition of grassland bacterial communities recognize the dominance of αand ß-classes of Proteobacteria, Acidobacteria, Actinobacteria, Verrucomicrobia, Bacteroidetes, Chloroflexi, Planctomycetes, and Firmicutes [1,2,[54][55][56][57][58]. This is in line with our findings. In addition, we identified sensitive bacterial taxa negatively responding to cattle overwintering: Acidobacteria, Planctomycetes, α-Proteobacteria, Verrucomicrobia, and Elusimicrobia. In contrast, Actinobacteria, Bacteroidetes, Deinococcus-Thermus, Firmicutes (both Bacilli and Clostridia), Chloroflexi, and OD1 were favoured in cattle impacted soils. In a previous study [59], the authors additionally observed an increase in the abundance of ß-Proteobacteria with increasing cattle impact, and a decrease in the abundance of Actinobacteria. The discrepancy in the trends observed may be associated with the different approaches used in the two studies, caused mainly by the PCR primer bias [60,61]. Bacterial diversity and community structure in soils is known to be mediated by vegetation and soil chemistry (mainly pH, organic matter, and nitrogen content). Bacteria, for example, have a narrower pH growth optimum than fungi [62], hence they are more sensitive to fluctuations in soil pH. An increase in the abundance of some bacterial groups in the soil, therefore, may reflect shifts in soil pH or nutrient content. O'Callaghan et al. [6] reported an increase in Firmicutes following bovine urine amendment. Nacke et al. [1] showed strong positive correlations between Bacteroidetes, Actinobacteria, and ß-Proteobacteria with soil pH, and a negative response of α-Proteobacteria. Indeed, differences in soil bacterial composition as a result of changes in soil chemistry caused by differing grazing intensity have been found by a number of authors [9,10,63,64]. However, most of them did not provide comprehensive data on the bacterial community as they reported community profiles based on 16S rRNA gene fingerprinting or PLFA. In our previous study [9], shifts in bacterial community caused by cattle use were recognized using DGGE community profiling, with cattle impacted soils clustered together apart from the control, which corroborates the results obtained by pyrosequencing (Fig 1B). A laboratory experiment carried out with arable soils showed shifts in bacterial community composition following pig manure amendment [58]. This experiment indicated that the application of pig manure had a rapid impact on the bacterial community, with a reduced abundance of Acidobacteria and Planctomycetes and an increase in Firmicutes after three days after pig manure addition. After 60 days, however, the effects of manure addition were marginal, thereby indicating a relatively fast recovery and high resilience of the bacterial community. In the field study outlined here, however, long-term effects of cattle overwintering on the bacterial community provoked changes similar to those immediately following manure addition in the laboratory [58], with similar trends in sensitive (Acidobacteria and Planctomycetes) and introduced (Firmicutes) bacterial groups. The study of Ding et al. [58], however, did not include the long-term effects of laboratory manure application alongside a field study. To date, no field study had been performed at sites with long-term and repetitive deposition of large amounts of manure during winter season (in our study, the LTI site receives manure for almost six months each year). Field studies of this type are important as land-use changes result in more complex effects than in laboratory studies, i.e. the combination of a range of environmental variables results in changes to vegetation and the topsoil (physical damage) and changes in soil physicochemical properties. The REG soil in our study showed significant differences in bacterial community composition compared to the LTI soil, indicated by the dominance of Acidobacteria and reduced anaerobic indicator groups (i.e. Clostridiales, Chloroflexi, and Bacteroidetes). In addition to changes in soil physicochemical properties, REG soil bacteria may also be strongly influenced by plant revitalization and increased soil aeration. Indeed, the connection between vegetation cover and soil bacteria community structure has been well studied [1,65,66]. Apparently, REG soil is in a transitory stage and its archaeal and bacterial communities more resemble those of CON than those of LTI.
On a global scale, Acidobacteria (an oligotrophic taxon) prefer low N content [67]. Similarly, Acidobacteria abundance decreased in cattle-impacted soils, with lowest levels observed in LTI soil. This quantitative change was also linked with qualitative changes. The highest Acidobacteria diversity was observed in CON soil, with the main groups represented being 1, 2, 3, 4, 5, 6, and 7. In cattle impacted soils, abundance of Groups 1, 2, 3, 5, and 7 decreased in response to altered soil conditions. Our findings could indicate that Groups 4 and 6 are more resilient to changes in pH, because soil pH ranged more broadly (5.2-7.9; Table 1) than that in the previously reported study (6.0-7.3; [2]). Though few bacterial indicators of cattle manure contamination were recognized at the genus level (S3 Table, Butyrivibrio, Clostridium, Peptostreptococcaceae Insertae Sedis, and Turicibacter) all increased in abundance in cattle impacted soils (mainly LTI and STI soil). Additionally, a number of other genera increased in cattleimpacted soils as a result of altered edaphic factors (Acidovorax, Caldilinea, Corynebacterium, Devosia, Dokdonella, Hirschia, Hydrogenophaga, Limnobacter, Nocardioides, Opitutus, Propionibacteriaceae bacterium, Proteiniclasticum, Pseudoxanthomonas, Smithella, Tetrasphaera, Thaurea, Trichococcus, and Truepera, S3 Table), which may indicate their high nutritional demands. Most of these bacteria are anaerobic or facultative anaerobic and the enrichment of Bacteroidetes (obligate anaerobes) in LTI soil in particular indicates an increase in the availability of anaerobic niches as a result of intensive and long-term cattle husbandry. Elhottová et al. [10] were able to show an increase in anaerobic microbial biomass in soils intensively managed winter pastures and, based on these and recent findings, we hypothesize that some rumen-borne microbes may survive long-term in such soils.
The significant increase of Chloroflexi in LTI soil is not in line with previous findings, which tend to show a preference for oligotrophic rather than nutrient rich soils [2,56,67]. This might be explained by changes in anaerobic status and/or through a higher availability of electron acceptors such as organohalogens [68]. Quantities of strictly and/or facultative anaerobic bacteria (e.g. Clostridia, Bacteroidetes, and Tenericutes) are transferred to the soil following deposition of manure, each of which have different abilities to survive in soil. Whereas Clostridia and Bacteroidetes clearly can survive in cattle-impacted soils ( Table 4), Tenericutes (Mollicutes) cannot as they are strict anaerobes preferentially associated with eukaryotic hosts and lack peptidoglycan cell wall.
Our results reflect the response of soil archaeal and bacterial communities to changes in cattle outdoor husbandry, further studies at similar overwintering farms in different localities are needed to generalize observed patterns and link them to changes in soil functioning.

Conclusions
Archaeal and bacterial community changes in upland grassland soils were related to different levels of cattle impact. The communities responded differently to environmental change, with archaeal diversity decreasing significantly in cattle-impacted soils and bacterial diversity less affected. Both prokaryote groups showed shifts in community composition, the shift being more pronounced in the Archaea. Methanogenic archaea and MCG increased in abundance in cattle-impacted soils at the expense of the soil indicator group Thaumarchaeota. Similarly, Clostridia, Actinobacteria, Bacteroidetes, and RF3, which were all present in cattle manure, were all enriched in cattle-impacted soils, while Acidobacteria were reduced. Additionally, several soil born archaea and bacteria (e.g. Methanosarcinaceae and Chloroflexi) multiplied under altered soil conditions. In general, anaerobic archaea and bacteria were identified as indicator groups for intensive and long-term cattle grazing. Of these, the methanogens and syntrophic bacteria indicated a preference for anaerobic processes, resulting in methane emissions from such cattle-affected soils. Recovery from cattle husbandry was characterized by a return to Acidobacterial dominance and the reduction of Clostridiales and other anaerobic groups, possibly related to revitalization of the vegetation cover.  Table. Mean relative abundance of bacterial genera in soils (LTI = long-term impact; REG = regenerating soil; STI = short-term impact; CON = control) and manure (CMN). Taxonomy was assigned using LCA Classifier against the SilvaMod database. Gradient of grey scale indicates the relative abundance (the darker color the higher abundance). (DOC) S4 Table. Results of Kruskal-Wallis ANOVA and Spearman's rank correlations between relative abundances of taxonomic groups (phyla, family or genera) and soil properties. All correlations shown are significant at P < 0.05; significance at P < 0.01 is indicated by underlining and bold type. Only the 20 most abundant bacterial genera in the dataset were tested. (DOCX)