Differences in microbial community structure and nitrogen cycling in natural and drained tropical peatland soils

Tropical peatlands, which play a crucial role in the maintenance of different ecosystem services, are increasingly drained for agriculture, forestry, peat extraction and human settlement purposes. The present study investigated the differences between natural and drained sites of a tropical peatland in the community structure of soil bacteria and archaea and their potential to perform nitrogen transformation processes. The results indicate significant dissimilarities in the structure of soil bacterial and archaeal communities as well as nirK, nirS, nosZ, nifH and archaeal amoA gene-possessing microbial communities. The reduced denitrification and N2-fixing potential was detected in the drained tropical peatland soil. In undisturbed peatland soil, the N2O emission was primarily related to nirS-type denitrifiers and dissimilatory nitrate reduction to ammonium, while the conversion of N2O to N2 was controlled by microbes possessing nosZ clade I genes. The denitrifying microbial community of the drained site differed significantly from the natural site community. The main reducers of N2O were microbes harbouring nosZ clade II genes in the drained site. Additionally, the importance of DNRA process as one of the controlling mechanisms of N2O fluxes in the natural peatlands of the tropics revealed from the results of the study.

Helium atmosphere soil incubation technique 3,4 was used to measure potential N 2 fluxes from soil cores in the same laboratory. The cylinders with the intact soil cores were placed into special gas-tight incubation vessels locating in the climate chamber. Gases were removed by flushing with an artificial gas mixture (21.0% O 2 , 358 ppm CO 2 , 0.313 ppm N 2 O, 1.67 ppm CH 4

Preparation of DNA libraries, sequencing and data processing
A quality check on raw sequence data was performed using FastQC v. 0.11.4 5 . About 1% of pairedend reads had low quality caused by sequencing errors. Quality trimming to remove poly-G tails and reads that have ambiguous nucleotides was done with Cutadapt v. 1.9.1 6 . DNA yield and characteristics of the sequencing data are shown in Supplementary Table 7. Kaiju v. 1.4.5, which search strategy finds maximum exact matches on the protein-level between query and a reference database (NCBI RefSeq database) using the Borrows-Wheeler transform, was used to classify metagenomic reads down to the species level 7 . 39.8±0.6% (natural site) and 34.3±3.2% (drained site) of the reads were classified with the Kaiju in the heuristic Greedy mode with the default parameter choices, allowing up to five amino acid substitutions during the search.
To screen the metagenomes of the samples for potential to perform nitrogen transformation processes, existing databases of marker genes (amino acid sequences) were used as a reference: nirK 8 ; nirS and nosZ 9 ; nifH, nrfA, hzsA, bacterial amoA/pmoA and archaeal amoA 10 . Additionally, phylogenetic marker rpoB gene (encoding the RNA polymerase) was used as a single copy gene reference 10 . The reference gene datasets alignments were generated using MUSCLE v. 3.8.31 11 and the respective phylogenetic trees were built for detected genes using FastTree v. 2.1.3 12 . Prodigal v. 2.6.2 was used to predict protein-coding regions 13 . Identification of the nitrogen cycling genes from metagenomes was performed with a trained HMM profile for each functional gene using hmmsearch from HMMER v.
ART v. GreatSmokyMountains-04-17-2016 was used to generate simulated sequencing reads dataset based on reference data in order to specify the statistical significance thresholds (E-values of hmmsearch) for limiting false positive matches against reference sequences for each nitrogen cycling marker gene (0.01 for nirK, nosZ and archaeal amoA; 0.001 for nirS, nifH, nrfA and hzsA; 0.0001 for bacterial amoA/pmoA) 17 . The sensitivity and specificity of the tests were on average 81.6% and 96.7%, respectively. References were tested against different homologs to control their possibility to pick up non-specific sequences and problem occurred in case of nirK and nirS references, where there is possibility that some of the hits may be homologs and these genera were marked on the phylogenetic trees of nirK ( Supplementary Fig. 1) and nirS ( Supplementary Fig. 2).
The overall percentages of the reads matching nitrogen cycle related genes were 28.2±2.6% (natural site) and 29.4±6.0% (drained site). Nitrogen cycling gene read abundances are shown as proportions of total rpoB (RNA polymerase) gene read abundance and normalised according to gene length ( Figure   4). Abundance-weighted phylogenetic diversity (balance-weighted phylogenetic diversity, BWPD) index values were calculated according to McCoy and Matsen 18 setting Ɵ=0.5. The variation in the composition of the N-cycling bacterial and archaeal communities among the study sites (beta diversity) was decomposed into replacement, richness difference, and nestedness components 19 .
Bacterial and archaeal phyla (profile obtained by metagenomics) abundances according to copy numbers per gram of dry soil of bacterial and archaeal 16S rRNA gene (obtained by qPCR) were shown in Supplementary Fig. 9.
Stock solutions of target sequence containing plasmids or in case of nifH gene sequence containing synthetic double-stranded DNA fragments (Eurofins MWG Operon, Ebersberg, Germany) were used to create serially diluted standard curves ranging from 25 to 10 9 copies for each target gene (Supplementary Table 9 The abundance of total prokaryotic organisms was calculated by summing the bacterial and archaeal 16S rRNA gene abundances. The proportion of target genes and microbial groups in soil microbial community was estimated by normalisations against prokaryotic community. The ratios between two types of nir genes (nirS/nirK), two nosZ gene clades (nosZI/nosZII) and nir and nosZ genes (nosZ/nir) were also calculated.

SUPPLEMENTARY TABLES
Supplementary