Short-term effect of Eucalyptus plantations on soil microbial communities and soil-atmosphere methane and nitrous oxide exchange

Soil greenhouse gas (GHG) emissions are a significant environmental problem resulting from microbially-mediated nitrogen (N) and carbon (C) cycling. This study aimed to investigate the impact of Eucalyptus plantations on the structure and function of a soil microbial community, and how resulting alterations may be linked to GHG fluxes. We sampled and monitored two adjacent Eucalyptus plantations—a recently logged site that harbored new seedlings and an adult plantation—and compared them to a site hosting native vegetation. We used 16S rRNA gene sequencing and qPCR amplifications of key nitrogen and methane cycle genes to characterize microbial structure and functional gene abundance and compared our data with soil parameters and GHG fluxes. Both microbial community attributes were significantly affected by land use and logging of Eucalyptus plantations. The genes nosZ and archaeal amoA were significantly more abundant in native forest than in either young or old Eucalyptus plantations. Statistical analyses suggest that land use type has a greater impact on microbial community structure and functional gene abundance than Eucalyptus rotation. There was no correlation between GHG fluxes and shifts in microbial community, suggesting that microbial community structure and functional gene abundance are not the main drivers of GHG fluxes in this system.

Of these Brazilian planted forests, 5.56 million ha are dedicated to Eucalyptus. Eucalyptus are fast-growing trees with high carbon sequestration potential during development 8,9 . Eucalyptus plantations have been reported to be a source of N 2 O and CO 2 and a sink of CH 4 in semi-arid and subtropical climates, as observed for most forest ecosystems 10,11 . However, GHG fluxes in Eucalyptus plantations have not yet been well described in the tropics, so a greater understanding of the impacts of Eucalyptus plantation management on these fluxes is still needed.
Soil behaves as both source and sink for GHGs 12 , as it represents the living space for the microbial communities responsible for nutrient cycling 13 . Accordingly, microbial activities in the N and C cycles are central to GHG fluxes in soil.
The link between soil microbial communities and GHG fluxes has previously been described 14 . Soil microbial processes are particularly impacted by land use practices, which can deregulate nutrient cycles and thereby increase or reduce GHG emissions 15,16 . Some studies have revealed a correlation between the abundance and/or expression of functional genes involved in N and C cycles and GHG fluxes in forest soils [16][17][18] . However, until now, no study has focused on the link between a microbial community and the GHG fluxes in the soil of Eucalyptus plantations.
To address this topic, we studied GHG fluxes and the microbial community associated with Eucalyptus plantations at two growth stages (i.e., one with new seedlings and one with 6-year-old trees), and with a native Brazilian tropical forest (Atlantic Forest). We hypothesized that: (1) replacement of native vegetation by Eucalyptus plantation or Eucalyptus plantation rotation would lead to changes in microbial community structure and functional gene abundance; and (2) changes in microbial community would alter GHG flux dynamics at each site.
We employed two different experimental strategies to test these hypotheses. Firstly, we performed 16S ribosomal RNA (rRNA) gene sequencing to compare the composition of soil microbial communities among the three treatments (i.e., young and old Eucalyptus plantations, as well as native forest). Secondly, we measured by quantitative PCR (qPCR) the abundances of the following key functional genes involved in the methane and N cycles: nifH (nitrogen fixation), archaeal and bacterial amoA (nitrification), nirK (denitrification), nosZ (N 2 O reduction) and mcrA (CH 4 production). We then compared our 16S rRNA sequencing and qPCR data with soil physicochemical properties and measurements of GHG fluxes.

Materials and Methods
Field site description. The experimental field site is located in Belo Oriente, Minas Gerais, Brazil (18-20°S, 42-44°W, 300 m elevation), in an area that belongs to the Celulose Nippo-Brasileira (CENIBRA) company. The climate is classified as Aw (tropical with a dry winter and a rainy season during summer, Köppen classification), with an annual mean temperature varying from 22 °C to 27 °C (maximum = 32 °C, minimum = 18 °C). Annual mean precipitation varies from 701 to 1,500 mm. The soil is classified as loamy red-yellow Ferralsol. The experimental area is highly sloping, with a slope of 26 degrees.
The experimental field was originally covered by Atlantic Forest, a native tropical forest. Since 1960, CENIBRA has managed Eucalyptus plantations in this area and adopts regular rotation cycles of 7 to 9 years between planting seedlings and tree-cutting. Part of the area has been retained as native vegetation, in accordance with Brazilian law, which allowed us to compare adjacent areas covered by native forest or Eucalyptus plantation. Planted seedlings are clones of Eucalyptus urograndis produced by the company.
We chose an area under Eucalyptus plantation since 1978, immediately adjacent to a fragment of Atlantic Forest to perform this study. The most recent Eucalyptus rotation started in 2011, with trees in the 6 th year growth in the beginning of 2017. To understand the short-term impact of a new rotation, we manually logged approximately half of this 6-year-old Eucalyptus plantation in February 2017 and replanted it immediately with new Eucalyptus seedlings. Sampling for our analyses was conducted at the end of March 2017. Experimental Design. Three adjacent areas under contrasting use were considered for this study: i. Atlantic Forest fragment (NF) -native Brazilian tropical forest fragment (Atlantic Forest); ii. Old Eucalyptus (OE) -Eucalyptus plantation with 6-year-old trees; iii. Young Eucalyptus (YE) -Eucalyptus plantation in which 6-year-old trees had been removed from the field, and new seedlings were planted one month before sampling.
We compared the Eucalyptus plantations at the beginning (YE) and at the end of the rotation cycle (OE), with the Atlantic Forest (NF) fragment acting as a reference. OE and YE areas comprised 470 and 680 Eucalyptus trees, respectively, which had been planted in lines with a spacing of 3 × 2.5 m. GHG sampling. Gas sampling was done daily, during four days, from 14 to 17 March 2017. Nitrous oxide and methane fluxes were manually quantified using closed static chambers, similar to those described by Alves et al. 19 . The top-base chambers had a base composed of a rectangular steel frame (40 cm × 60 cm), which was deployed between rows of Eucalyptus trees or randomly in the Atlantic Forest fragment. The base was inserted into the soil to a depth of 6-7 cm, before attaching a polyethylene lid (of the same lateral dimensions as the base) to it, generating an internal chamber space 12-15 cm above the soil surface. Soft rubber was fixed to the rim of the lid to satisfactorily seal it to the base. Air accumulating in the chamber headspace was withdrawn through a three-way valve connected to the lid using a polypropylene syringe. Approximately 30 mL of this air was sampled and transferred to 20 mL chromatography vials crimped with chlorobutyl septa. Just before gas transfer, each vial was evacuated to approximately −100 kPa by using an electric vacuum pumping system. The chamber lid was covered with a 2-cm-thick foam layer and reflective adherent mantle for thermic insulation. The bases of the chambers remained in their respective areas throughout the monitoring period. After chamber closure an air ScIEntIFIc RepoRTS | (2018) 8:15133 | DOI:10.1038/s41598-018-33594-6 sample was immediately taken from chamber headspace. Subsequent samples were taken at 20, 40 and 60 min, after which lids were removed. Gas sampling was always performed in the morning between 08 h and 10 h 19 .
Gas analyses were performed using a gas chromatograph (GC 2014, Shimadzu, Japan). For each sampling run, N 2 O and CH 4 standards were used to build an analytical curve to transform the integrated areas of each sample peak into gas concentrations.
Gas fluxes were calculated by the equation F = (ΔC/Δt)(V/A)M/Vm, where ΔC/Δt is the slope of a linear function fitted to the gas concentration of samples taken at 0, 20, 40 and 60 min after chamber closure; V and A are the volume of the chamber and the area of soil covered by the chamber, respectively; M is the molecular weight of atoms of the elements N and C, respectively, in molecules of N 2 O and CH 4 ; and Vm is the molecular volume at the sampling temperature.
Biogeochemical characteristics of the soil. We sampled approximately 500 g of soil from nearby each of the gas chambers described in the previous section (around 10 cm from the chamber) to measure clay content, pH, total carbon, total nitrogen, available P, and amounts of H + + Al 3+ , Ca 2+ , Mg 2+ and Al 3+ . Each sampling point served as an independent replicate, resulting in five replicates per treatment, which were used to assess correlations with the respective gas fluxes calculated for each chamber.
Soil samples were air dried, sieved (2 mm) and analysed chemically. Total carbon and total nitrogen were determined using a CHN elemental analyser (PerkinElmer, USA). Exchangeable nutrients: Ca 2+ , Mg 2+ and Al 3+ were extracted by 1 M KCl; P, Na and K by Mehlich-1 extractant −0.05 mol L −1 in HCl in 0,0125 mol L −1 H 2 SO 4 ) and pH (soil:water, 1:10); Potential acidity: H + Al extracted with calcium acetate 1 N (pH 7), titrated with 0.0125 NaOH N. Inductively coupled plasma apparatus for Ca +2 , Mg +2 and Al + , flame emission (K and Na) and photocolorimeter (for P) were used to nutrient determinations. Soil granulometry was determined using the aerometer method, after chemical dispersion. All these soil characteristics were measured according to Embrapa 20 .
Soil moisture and concentrations of NO 3 − and NH 4 + were measured according to Morais et al. 21 . Determinations were made from samples taken from another four randomly selected points in each of the three areas. The inorganic N content was extracted from fresh soil with 60 mL 2M KCl after 1 h on a rotary shaker at 220 rpm. The supernatant was filtered and the NO 3 − and NH 4 + concentrations were determined in the resultant solution respectively by flow injection (FIA) technique using Cd reduction and nitrite analysis and by the salicylate reaction adapted for FIA. We used averages of the four values obtained for soil moisture and NO 3 − and NH 4 + concentrations.
Soil sampling for microbial analyses. Soil samples were taken from five points in each treatment area (YE, OE and NF), approximately 10 cm away from each gas flux chamber, resulting in five replicates per treatment. To extract each soil sample, we placed a steel tube probe of 1.5 cm diameter (previously sterilized at 180 °C for 3 h to remove contaminants and nucleases) into the soil to a depth of 7 cm. The harvested soil was immediately put in a sterile 50 mL propylene tube for mixing, before being separated into two subsamples and placed in liquid nitrogen until we conducted DNA extractions.
Analysis of bacterial community structure. DNA was extracted from approximately 500 mg of each soil sample using the Fast DNA Spin Kit for Soil (MP Biomedicals, USA). The extracted DNA was dark in color due to a high level of humic material, so it was then purified using the last steps of the NucleoSpin Soil kit (Macherey-Nagel, Germany) protocol. DNA concentration and quality were measured using a NanoDrop (Thermo Fisher, USA). Soil bacterial diversity were assessed by next generation sequencing of the V4 variable region of the 16S rRNA gene using the primers 515FB (GTGYCAGCMGCCGCGGTAA) and 806RB (GGACTACNVGGGTWTCTAAT) 22,23 . PCR reactions with a barcode on the forward primer were used in a 28-cycle PCR using the HotStarTaq Plus Master Mix Kit (Qiagen, USA) under the following conditions: 94 °C 3 min; 28 cycles of (94 °C 30 s, 53 °C 40 s, 72 °C 1 min); 72 °C 5 min. Following PCR amplification, PCR products were checked in 2% agarose gels to determine the success of PCR and the relative intensity of resulting bands. Multiple samples were pooled in equal proportions based on their molecular weight and DNA concentrations. Pooled samples were purified using calibrated Ampure XP beads (Beckman Coulter, USA). The pooled and purified PCR products were then used to prepare an Illumina DNA library.
Sequencing was performed by MR DNA (www.mrdnalab.com, Shallowater, TX, USA) on a MiSeq (Illumina, USA) paired-end 2 × 250 sequencing system, following the manufacturer's guidelines. Raw data were downloaded from Basespace and analysed using Mothur v1.39 software 24 . Forward and reverse paired sequences were merged into contigs after checking for the presence of barcodes and primers in the sequences. Merged sequences of less than 240 base pairs (bp) or greater than 300 bp, containing any ambiguities, or containing more than 8-mer homopolymers were removed.
Sequences were then aligned using a modified Silva database (generated by a virtual PCR using the same primers as those used for our samples) as reference 25 , and the resulting alignment was submitted to screen. seqs and filter.seqs (Mothur v1.39) to remove badly aligned sequences and uninformative columns in the alignment. The sequences where then pre-clustered using the command pre.cluster (Mothur v1.39.5) with parameter diffs = 2. Chimeras were detected using the command chimera.vsearch (Mothur v1.39.5) and then eliminated. We classified sequences using the classify.seqs (Mothur v1.39.5) command, with the Ribosomal Database Project 26 as reference and a bootstrap threshold of 80. Sequences identified as being from chloroplasts, mitochondria, Eukarya or Archaea and those not assigned to any kingdom were removed. The resulting sequences were used as input for the dist.seqs (Mothur v1.39.5) command. Finally, we clustered the sequences into operational taxonomic units (OTUs), with a 3% dissimilarity threshold. To avoid bias due to sampling effort, the samples were randomly normalized to the same number of sequences (35028). We employed a taxonomic summary to assess the bacterial composition of each sample.
Differences in the relative frequencies of phyla and classes among treatments were tested using analysis of variance (ANOVA) followed by Tukey's test. Before performing ANOVA, we checked the homoscedasticity among treatments and normality of distributions (Shapiro-Wilk's test) for all values, and data was transformed (log or Box-Cox transformation) accordingly if necessary before performing ANOVA. The Shannon index was also analysed by ANOVA followed by a pairwise Tukey's test.
The distribution of OTUs was used as input for a non-metric multidimensional scaling (NMDS) ordination with the Bray-Curtis similarity index to assess relationships among samples. We performed a PERMANOVA test, followed by Bonferroni correction (p < 0.05), to assess differences in structural composition among treatments. All aforementioned statistical analyses were conducted in Past3 27 . To determine if the treatments had a significant effect on specific bacterial OTUs, we undertook a blocked Indicator Species Analysis (ISA) (YE vs OE and YE + OE vs NF) using Mothur v1.39.5. The ISA, proposed by Dufrêne and Legendre 28 , is based on the relative frequency of a specie (OTU in our case) within and inter treatments. It gives provides an indicator value ranging from 100 (perfect indicator) to 0 (no indicator). The perfect indicator species would be the one found in high abundance in all samples of a given treatment and absent in all other treatments. It also uses a randomization test to evaluate the significance of the specie distribution.
Nucleotide sequence accession numbers. The data generated were deposited in the NCBI Sequence Read Archive (SRA) and are available under Bioproject accession numbers PRJNA471919.
Gene quantification by qPCR. Genes were chosen based on their involvement in the soil nitrogen and methane cycles. Primers were selected according to the literature and assessed with the Primer blast tool of the National Center of Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). We chose primers with the highest number of results and lacking non-specificity ( Table 1).
The standards were constructed by amplifying each gene from DNA extracted from soil or activated sludge, ligating them into plasmids (CloneJET PCR Cloning Kit, Thermo Fisher Scientific, USA), and transforming them into E. coli DH5alpha plasmid (heat shock method, Froger and Hall, 2007). The plasmids were recovered using the PureYield Plasmid Miniprep system (Promega, USA). Based on the size of each gene, the weight of one nucleotide, and the plasmid concentration, we generated ten-fold serial dilutions in RNase-and DNase-free water for each plasmid to reach 10 10 to 10 2 gene copies per reaction.
qPCR reactions were carried out using the GoTaq ® qPCR Master Mix (Promega, USA) on DNA (DNA amount ranged from 56 to 106 ng) in a 7500 Real-Time PCR system (Applied Biosystem, USA) with the SybrGreen excitation setting.
Each sample was quantified twice (technical replicate) in a 20 μl reaction volume. The following program was used: 50 °C 2 min; 95 °C 5 min; 40 cycles of (95 °C 30 s, annealing temperature (Table 1)  cycle. An absolute quantification was realized based on the standard curve generated by the plasmid dilutions. The copy number of each sample was normalized to the weight of soil used for DNA extraction to take into account the difference in soil DNA abundance among treatments. qPCR efficiencies were calculated using the formula E = 10 (−1/slope) − 1, (where slope is the slope of the standard curve), with E = 1 corresponding to 100% efficiency.
Each qPCR reaction result was subjected to ANOVA followed by a pairwise Tukey's test. We ran an NMDS ordination based on soil characteristics, using qPCR reaction results and gas fluxes as correlating parameters. Spearman correlations were generated between gene amounts, gas fluxes and soil characteristics. We only present significant correlations (p < 0.05).

Results
Soil characteristics. Our soil texture analysis revealed similar clay contents among all three studied areas (Table 2). However, we did identify some differences in soil characteristics among treatments, especially lower soil pH in area YE than NF, and higher NO 3 − content in areas YE than OE compared to NF. Area YE exhibited higher intra-sample heterogeneity than OE and NF. Despite high spatial variability, the low levels of available phosphorus in soil under native vegetation (area NF) reveal the importance of fertilization to achieve better wood production in high weathering soils. Fertilization resulted in considerable differences in the N:P soil ratio.
Gas fluxes. All three treatment areas exhibited similar GHG flux dynamics (Table 3) over the short experimental period. All three areas acted as a sink for CH 4 , with average fluxes ranging from −22 to −35 μg m −2 h −1 . All three areas acted as sources of N 2 O, with average fluxes ranging from 4.8 to 9.4 μg m −2 h −1 . These fluxes did not differ significantly among treatments.

Impact of land use change and Eucalyptus logging on microbial community structure.
We assessed a total of 525,420 sequences (35,028 per sample) after applying quality filters and data normalization, clustered into 6,831 OTUs (3% dissimilarity threshold). Rarefaction curves show that our sequencing effort describes well the diversity of each sample (Supplemental Fig. 1). Bacterial richness (represented by numbers of OTUs) was significantly different among treatments, being highest in YE, followed by OE, and lastly NF (Table 4). Shannon diversity indices also differed significantly among treatments, being higher in YE and OE than in NF (Table 4). No significant difference was found between the two Eucalyptus plantations.
Taxonomic assignments revealed that all treatment areas were dominated by the same phyla, but with significant differences in the abundances of some phyla among treatments (Fig. 1A). Eight phyla were found in all treatments: Proteobacteria, Acidobacteria, Actinobacteria, Verrucomicrobia, Chloroflexi, Firmicutes and Bacteroidetes. The three most dominant phyla (Proteobacteria, Acidobacteria and Actinobacteria) represented at least 74% of the communities in each of the three treatments. We could classify approximately 88% of sequences into 16 different classes (Fig. 1B). Alphaproteobacteria were dominant in all treatments, representing 25 to 32% of the entire bacterial community. Most significant differences in the relative abundances of bacterial taxa were observed between the native forest treatment (NF) and the two Eucalyptus plantations (YE and OE). However, there were also significant differences in the microbial community between the YE and OE treatments.
A global analysis of the distribution of the 6,831 OTUs indicated that, overall, microbial community structure was significantly affected by treatments, with all treatments differing from one another (PERMANOVA, p < 0.001). To identify the main bacterial groups responsible for the structure change, the abundances of the   50 most abundant OTUs (representing 56% of the entire community) were tested by a ISA, which allowed us to identify 23 OTUs that were significantly different between the YE + OE and NF treatments. The 23 OTUs identification and their relative abundances in each treatment were shown in Fig. 2A. Additionally, the ISA performed on the 50 most abundant OTUs present in both Eucalyptus treatments revealed only six OTUs that were significantly different between them. The six OTUs identification and their relative abundances are shown in the Fig. 2B.

Quantification of abundances of C and N cycle functional genes.
Our qPCRs were efficient in terms of quantifying functional gene copy numbers (Fig. 3). Dissociation curves indicated that the reactions were specific for all genes (data not shown), with R-squared values of the standard curves ranging from 0.98 to 0.99. Copy numbers of the 16S rRNA gene were very similar among treatments (approx. 10 10 copies). No significant differences were observed among treatments for the genes nifH (approx. 10 9 copies), nirK (approx. 10 7 copies), bacterial amoA (approx. 10 9 copies), or mcrA (approx. 10 7 copies). However, copy numbers of the genes nosZ and archaeal amoA were significantly affected by land use, with both being lower in the two Eucalyptus plantations than in native forest. There were no significant differences between YE and OE for these two genes. Archaeal  amoA, bacterial amoA, nifH and nosZ (approx. 10 9 to 10 10 copies) were more abundant compared to nirK and mcrA (approx. 10 7 to 10 8 copies).
Correlations between abundances of N cycle genes, microbial community composition, soil characteristics, and gas fluxes. To more clearly understand how microbial activities are related to soil physicochemical characteristics, we generated correlations between abundances of N cycle genes (quantified by qPCR) and soil characteristics using Spearman correlations. Three genes exhibited correlations with soil characteristics-nifH, archaeal amoA and nosZ ( Table 5)-all of which were correlated with the N:P ratio and available P, whereas archaeal amoA and nosZ were both also correlated with humidity, NH 4 + and the C:N ratio. We also employed Spearman correlations to investigate the factors linked with gas fluxes in the three treatments, but there were no correlations between N 2 O and CH4 fluxes and N cycle gene abundances or any soil characteristic.
We ran non-metric ordinations (NMDS) to better visualize the data structure and differential relationships in terms of microbial community composition (OTUs) and functional gene abundance (qPCR). The ordination based on soil characteristics revealed that all treatments differed from each other and that the difference between the YE and OE treatments was more pronounced than the difference between the OE and NF treatments (Fig. 4A). This ordination also revealed higher variability within the YE treatment. An ordination based on the distribution of OTUs (Fig. 4B) indicated that differences between the NF and YE + OE treatments were stronger than those between YE and OE. These differences were correlated with gene copy numbers of archaeal amoA and nosZ (both of which were higher in the NF treatment) and with gene copy number of nirK, total C content and the C:N ratio (higher in the two Eucalyptus treatments). Moreover, the higher graphical dispersion for the YE and OE treatments compared to the NF treatment on the NMDS ordination indicates higher beta diversity in the Eucalyptus treatments than in the NF treatment.

Discussion
The main phyla we found in the soil samples of the three treatments are common in soils, either from cultivated or native forests [29][30][31][32][33] . A dominance of Proteobacteria has also previously been reported for many types of soils 30,31,33 . Despite the apparent homogeneity in microbial composition across treatments, the PERMANOVA on OTUs indicated that microbial community structure is significantly affected by both land use change (from native forest to Eucalyptus) and at the start of a new Eucalyptus rotation (the transition from the OE to YE treatment). Conversion of native forest to other land uses-such as silviculture, agriculture or pasture-has previously   been shown to have an impact on global soil microbial community structure 34,35 , but we are the first to report a short-term impact (within one month) of Eucalyptus plantation logging.
Our results suggest that the main driver of the observed differences in microbial community structure is land use type, as opposed to management practices (i.e., plantation logging and initiation of new rotations) or soil characteristics. Indeed, our NMDS analyses showed that OE and NF treatments were more similar to each other in terms of soil characteristics (Fig. 4A) relative to the YE treatment, but that the microbial community structure of OE was more similar to that of the YE treatment than the NF treatment, meaning that the overlying plant coverage was more effective than soil characteristics in driving microbial community structure (Fig. 4B). A similar outcome has been reported in other studies [36][37][38] , which indicate that this scenario is more likely to occur when soil pH varies little between treatments, as found in this study. Moreover, our comparison of relative abundances of specific OTUs between the YE + OE treatments and the NF treatment ( Fig. 2A) evidence that Eucalyptus plantations considerably influence the soil microbial community, including for specific groups of bacteria. High level selection of specific fungal groups by Eucalyptus trees has previously been reported 39 .
The higher bacterial alpha diversity (Table 4) and beta diversity (Fig. 4B) of the two Eucalyptus treatments compared to the native forest is unexpected considering the greater tree biodiversity of the latter treatment. Interestingly, increased alpha diversity was also reported for an Amazon forest zone 4 months after deforestation 35 . We postulate that microbial richness might increase under these circumstances to adapt to the disturbance caused by deforestation, which is supported by the higher diversity of the YE area in our study given that its tree coverage was removed one month before soil sampling.
However, the higher microbial diversity in the OE area cannot be linked to short-term disturbance, since the trees had been growing there for 6 years before sampling. Higher primary productivity in Eucalyptus plantations compared to native forest might explain the higher microbial diversity under Eucalyptus, perhaps leading to higher fluxes of root exudates into the soil. The higher soil beta diversity under Eucalyptus plantations relative to native forest might be explained by heterogeneous perturbations due to management practices, such as fertilization and the physical consequences of tree-cutting and -dragging before starting a new rotation.
The phyla we identified as being impacted by the three treatments (Proteobacteria, Acidobacteria, Actinobacteria, Planctomycetes, Verrucomicrobia; Fig. 1A) have previously been highlighted as being affected by altered land use, i.e., from forest cover to deforested 34,35,40 . The decreased relative abundance of Proteobacteria we observed between native forest and Eucalyptus plantations is consistent with other studies. A similar outcome was observed 20 to 30 years after native forest had been converted to oil palm plantation 34 , and also 2 to 3 years after an Amazonian forest soil lost its forest coverage 40 . Interestingly, both Acidobacteria and Actinobacteria were affected by starting a new Eucalyptus rotation (YE vs OE) but not by altered land use (NF vs YE + OE). This result suggests that plantation management practices may have their own inherent impact on microbial community structure.
We found the relative abundance of Acidobacteria to be significantly reduced in the YE treatment compared to that of the OE area. This phylum has previously been shown to be affected by deforestation 40 . Although relative abundances of Acidobacteria (particularly subdivision 1) are considered to be mainly driven by pH 41 , our results suggest that other factors can influence Acidobacteria abundance, as we did not find pH to be significantly correlated with Acidobacteria abundance (data not shown).
The homogenous distribution of 16S rRNA gene sequence abundances among treatments demonstrates that the size of the prokaryote population is not significantly different among treatments. Moreover, the abundance of 16S rRNA gene fragments in our study is within the range usually measured in soil 42 . The higher archaeal amoA gene copy number compared to bacterial amoA could be linked to the acidic soil pH (<4.5 in all treatments) 43 . This phenomenon has been widely observed in previous soil studies, and may be explained by the fact that NH 3 content is reduced at lower pH and ammonia-oxidizing Archaea have a higher affinity for NH 3 than ammonia-oxidizing bacteria 18 .
Since significant differences in gene abundances were only observed between native forest and Eucalyptus treatments, and not between the Eucalyptus treatments themselves, it seems that land use change and not starting a new rotation has an impact on microbial processes linked to the N and C cycles. Only the archaeal amoA and nosZ genes were significantly more abundant in soils under native forest than those under Eucalyptus. These differences are likely due to the specific soil characteristics of each treatment regime, since both archaeal amoA and nosZ gene abundances were strongly correlated with certain soil characteristics ( Table 5).
The negative correlation between archaeal amoA abundance and humidity might be due to the aerobic character of ammonia-oxidizing Archaea, even though these bacteria have been previously shown to tolerate a broad range of oxygen concentrations 44 . The negative correlation between the C:N ratio and archaeal amoA abundance seems logical, given that a lower C:N ratio means higher N abundance, so there is more substrate for ammonia-oxidizing Archaea 45 . However, a negative correlation was observed between NH 4 + and archaeal amoA abundance, suggesting that higher NH 4 + levels are more a consequence of low nitrifier activity than a cause of higher rates of nitrification 46 .
A similar process has been suggested to explain a negative correlation between soil NO 3 − levels and nirK gene abundance 38 . However, the correlations we report between levels of NH 4 + and NO 3 − and humidity should be viewed with caution, as these factors were not measured specifically at each soil-sampling site, but at four sites randomly distributed in the three treatment areas. Thus, the correlations are based on average values for each treatment area, which could lead to spurious correlations.
The negative correlation between nosZ gene copy number and humidity is unexpected, because nitrous oxide oxidase has been reported to be highly sensitive to O 2 , even being inhibited by low oxygen levels 47 . However, another study found a higher gene copy number of nosZ in wetlands than dry areas 48 , and a negative correlation between levels of nosZ and NH 4 + has been observed in wetlands 48 . This correlation might mean that lower NH 4 + levels are a consequence of higher NH 4 + oxidation activity from high levels of nitrous oxide reductase. It is important to note that we measured gene copy number, and not gene expression, in this study. Gene copy number is not always an effective means of reporting microbial activity as it provides no information about gene expression or protein activity.
The negative correlation we revealed between available P and the nifH, archaeal amoA and nosZ genes suggests an important link between P availability and N cycle functioning. However, this finding is not consistent with previous data, which suggest increased P availability promotes nitrification, denitrification and nitrogen fixation 49 . Analysing the bacterial community structure (Fig. 4) we can speculate that this correlation is more likely to be occurring indirectly. The microbial community structure of YE (highly disturbed) is positively correlated with P availability and negatively correlated with nifH and nosZ and AOA. Therefore, the disturbance could be the main factor controlling these correlations.
The absence of a correlation between levels of functional genes and gas fluxes is surprising because correlations between the nosZ, archaeal amoA, and bacterial amoA genes and N 2 O fluxes have previously been described in forest systems [16][17][18] . The quantification results for pmoA may explain CH 4 fluxes, as this gene has been shown to be positively correlated to CH 4 emissions from soil 17 . Moreover, we only conducted gene quantification for a single time-period, so we cannot rule out a correlation between gene abundances and gas fluxes over time. Despite the absence of an overall link between microbial structure and functional gene abundance and gas fluxes, it is important to note that our study focused solely on bacterial activities (apart from the archaeal amoA gene). It would be interesting to also analyse Archaea and fungi, both of which play important roles in N cycling in soils, especially in forest soils 50 .
In conclusion, we demonstrate that soils under Eucalyptus plantations can harbour a distinct and more diverse bacterial community compared to soils under native forest and that this community responds very quickly to environmental disturbances, such as implementation of a new plantation rotation. However, the short-term changes we observed arising from plantation management practices were qualitative and not quantitative, and they did not result in significant changes in terms of greenhouse gas fluxes from soil.